Active learning-enhanced neuroevolution potential for predictive modeling of UO₂ thermophysical properties

Abstract

Accurately characterizing the temperature dependence of UO₂ thermal conductivity is crucial for evaluating its performance under nuclear reactor operating conditions. However, experimental measurements are costly, density functional theory (DFT) calculations are constrained by small spatiotemporal scales, and traditional empirical potentials struggle to capture strong anharmonic effects. To this end, we developed a machine-learned neuroevolution potential (NEP) with near-DFT accuracy using an active learning strategy, and we systematically evaluated and cross-validated the thermal conductivity of UO₂ using equilibrium molecular dynamics (EMD), homogeneous nonequilibrium molecular dynamics (HNEMD), and nonequilibrium molecular dynamics (NEMD). The results demonstrate that HNEMD delivers a high signal-to-noise ratio, low uncertainty, and rapid convergence, exhibiting superior computational efficiency and robustness. At 500 K, the spectral phonon mean free path spans approximately one order of magnitude, and heat-transport channel lengths exceeding about 5 µm approach the bulk thermal conductivity limit. In the 800–1500 K range, the NEP reproduces the experimental temperature dependence of UO₂ thermal conductivity, while at lower temperatures (300–800 K), it achieves predictive accuracy comparable to that of DFT+U. Systematic validation of UO₂ fundamental properties including the equation of state, phonon dispersion relations, elastic constants, heat capacity, and linear thermal expansion coefficient demonstrates that the constructed NEP is reliable and broadly applicable. This work provides methodological support for multiscale thermal transport modeling of nuclear fuels and reactor safety assessment.

---

Introduction

Uranium dioxide (UO₂) is widely used as a nuclear fuel in commercial reactors due to its excellent thermal stability and high melting point. As a ceramic fuel, UO₂ exhibits heat transfer that is predominantly governed by quantized lattice vibrations (phonons).¹ The propagation of phonons within the material determines its thermal conductivity, which significantly influences the temperature distribution and thermal stresses within UO₂ fuel pellets.² These factors are critical to the safe and efficient operation of nuclear reactors. Therefore, accurate prediction of the thermophysical properties of UO₂, especially its thermal conductivity, is of paramount importance.

Traditionally, the thermophysical properties of UO₂ have been investigated using a combination of experimental and theoretical approaches. Numerous experimental correlations for predicting thermal conductivity have been proposed.^3–5 Fink et al.^6,7 systematically compiled experimental data on enthalpy, heat capacity, thermal expansion coefficients, and other properties. Their compilations remain highly valuable references to this day. However, experiments under harsh conditions such as high temperatures and intense radiation are technically demanding and costly. In addition, although theoretical approaches such as density functional theory (DFT) have been extensively applied to UO₂,^8–14 their computational expense limits simulations of large systems and long timescales. To address these limitations, empirical interatomic potentials (EIPs) have been developed.¹⁵ These potentials are typically built from physically analytical functions with parameters fitted to experimental data. Widely recognized EIPs for UO₂ include the Basak,¹⁶ Morelon,¹⁷ Yakub,¹⁸ and CRG¹⁹ potentials, among others. Using these potentials, one can study the properties of UO₂ over a wide temperature range such as lattice constants, bulk moduli, defect formation energies, and defect migration barriers. Nevertheless, due to their limited number of parameters and reliance on simplified analytical forms, EIPs often struggle to accurately capture complex interatomic interactions.

In recent years, machine learning interatomic potential (MLIP) has emerged as a promising computational tool. Trained on high-precision datasets, MLIPs can achieve near first-principles accuracy while significantly reducing computational costs. For instance, Dubois et al.²⁰ developed two types of MLIPs, namely the spectral neighbor analysis potential (SNAP) and a high-dimensional neural network potential (HDNNP), for simulating elastic properties, phonon dispersion relations, and defect formation energies in UO₂. Stippell et al.²¹ combined active learning with transfer learning to construct an MLIP for UO₂ that achieves DFT+U accuracy, enabling efficient prediction of properties such as lattice constants, bulk modulus, and radial distribution functions. Konashi et al.²² developed a neural network potential (NNP) for UO₂ and computed various thermophysical properties from room temperature to the melting point, including the molten state.

However, studies on MLIPs for UO₂ remain relatively limited, particularly in terms of systematic evaluations and public comparisons of thermal conductivity predictions. Moreover, most empirical potentials for UO₂ (e.g., Buckingham/Coulomb-types, CRG, etc.) were constructed without explicit constraints on phonon dispersion, Grüneisen parameters, or higher-order anharmonic interactions (such as three-phonon and four-phonon scattering). This often results in the predicted temperature dependence of thermal conductivity deviating from experiments and carrying substantial uncertainty.^23–25 Furthermore, the implementation of microscopic heat flux in molecular dynamics packages is method-dependent and may introduce biases unrelated to the potentials. For example, in the widely used large-scale atomic/molecular massively parallel simulator (LAMMPS),²⁶ the heat flux for many-body potentials typically relies on specific decompositions of atomic energy and stress, whose non-uniqueness generally precludes equivalence with the rigorous Hardy heat flux.^27–29 Within the Green–Kubo framework, using inconsistent heat-flux expressions can induce significant shifts in the computed thermal conductivity. By contrast, the graphics processing units molecular dynamics (GPUMD) package³⁰ developed by Fan et al. implements a Hardy-type heat flux that is consistent with energy conservation for many-body potentials.^31–33 It provides modules for thermal conductivity calculations via equilibrium molecular dynamics (EMD), homogeneous nonequilibrium molecular dynamics (HNEMD), and nonequilibrium molecular dynamics (NEMD). In addition, the neuroevolution potential (NEP) proposed by Fan et al.^34–36 has attracted widespread attention due to its exceptional computational efficiency. Nowadays, NEP has been successfully applied to a variety of material systems, including bulk crystals, two-dimensional materials, and amorphous phases, and has shown outstanding performance in thermal transport studies by combining GPUMD.

Therefore, we first generate an initial seed dataset of representative configurations and enrich it through active learning iterations to achieve sufficient sampling of the structural phase space across a wide temperature range, yielding a high-fidelity training dataset. On this basis, we train a NEP machine-learning potential. The reliability of the NEP is systematically assessed by comparisons with DFT for the equations of state, elastic constants, and phonon dispersion relations. Using the GPUMD package, we investigate the temperature dependence of the thermophysical properties of UO₂, including heat capacity and the linear thermal expansion coefficient. We provide accurate predictions of the thermal conductivity of UO₂ by cross-validating results from EMD, NEMD, and HNEMD, and we quantify the consistency among these methods.

Computational details

Construction of the initial dataset

To construct an initial training dataset for active learning, we adopted a stratified sampling strategy. First, using the NepTrainKit³⁷ software, we imposed lattice deformations ranging from −6% to 6% (including uniaxial, biaxial, and triaxial strains) on the UO₂ unit cell and a 2 × 2 × 2 supercell, together with a random atomic displacement of 0.1 Å. In total, 400 initial configurations were generated (313 UO₂ unit cells and 87 supercells). Subsequently, we augmented the dataset by applying uniform lattice expansions with scaling factors of 1.07–1.12 to the 2 × 2 × 2 supercell, with atomic displacement perturbations of 0.05 Å. Finally, using the NEP89³⁸ universal machine learning potential (covering 89 elements), we performed temperature-ramp molecular dynamics (MD) simulations on 3 × 2 × 2 and 3 × 3 × 2 supercells, heating from 1100 to 2500 K with a time step of 0.5 fs and a total duration of 7.5 ns. Representative configurations were then selected using the farthest point sampling (FPS) method integrated into NepTrainKit.³⁷ This approach systematically captures structural features from static deformations to dynamical evolution, ensuring dataset diversity while reducing redundancy.

Active learning

The configuration-space diversity of the training dataset and the accuracy of first-principles calculations are key determinants of the predictive accuracy of machine learning models. To enhance generalization capability and learning efficiency, we adopted an active learning (AL) strategy that iteratively selected and labeled high-value configurations to progressively improve model performance.

In this work, the active learning strategy consists of three steps: training, sampling, and labeling. (I) Training: four NEP models were trained with different random seeds on the initial dataset to enable subsequent configuration sampling and uncertainty estimation. (II) Sampling: one of these NEP models was randomly selected to perform MD simulations over a wide temperature range, generating diverse configuration samples. These configurations were then evaluated by all four models and those with high uncertainty indicated by disagreement among the model predictions were selected as candidates for labeling. (III) Labeling: the labeling criterion for candidate configurations is based on the maximum standard deviation of atomic forces,³⁰ denoted by σ_f:

where σ²_i,k is the sample variance of the predicted force for atom i along the Cartesian direction k, with k ∈ {x, y, z}. The labeling criteria were as follows: if σ_f ≤ 0.06 eV Å⁻¹, the configuration was considered sufficiently accurate and no additional DFT calculation was performed; if σ_f > 0.06 eV Å⁻¹, the configuration was flagged as high-uncertainty and DFT calculation was carried out to obtain reference labels. In addition, to ensure physical plausibility, we used the unphysical-structure identification tool in NepTrainKit³⁷ to remove highly distorted or energetically anomalous configurations (e.g., unrealistically short interatomic distances or unreasonable bond angles). High-accuracy single-point DFT calculations were then performed to obtain energies, atomic forces, and the virial tensors for the filtered structures, which were added to the training dataset. The NEP model was retrained on the updated dataset and iteratively refined to improve predictive accuracy and robustness. Active learning was terminated once fewer than 5% of newly sampled configurations exhibited σ_f > 0.06 eV Å⁻¹ and the predicted thermophysical properties stabilized over the target temperature range.

Across 13 rounds of active learning, we progressively optimized the quality and size of the training dataset by increasing MD simulation durations to more thoroughly explore the diversity and dynamical evolution of the configuration space. Ultimately, we constructed a high-quality dataset comprising 1918 independent configurations and 148 632 atoms, including 506 UO₂ unit cells, 1129 2 × 2 × 2, 162 3 × 2 × 2, 30 3 × 3 × 2, and 91 4 × 1 × 1 supercells. The settings for individual active learning iterations are summarized in SI, Table S1.

DFT calculation settings

All first-principles calculations were performed using the Vienna ab initio simulation Package (VASP).³⁹ In line with our focus on lattice thermophysical properties, which are dominated by lattice dynamics with negligible electronic contributions, and to ensure structural diversity and stable labels in active-learning, we used non-spin-polarized DFT without Hubbard corrections. Electron–ion interactions were described by the projector augmented-wave (PAW)⁴⁰ method, and exchange–correlation effects were treated within the generalized gradient approximation (GGA) with the Perdew–Burke–Ernzerhof (PBE)^41,42 functional. The 6s²7s²6p⁶6d²5f² and 2s²2p⁴ orbitals were treated as valence electrons for U and O, respectively. A plane-wave cutoff of 600 eV was employed, with Gaussian smearing (σ = 0.02 eV). The Brillouin zone was sampled on Γ-centered Monkhorst–Pack (MP)⁴³ grids. To balance sampling with the size of cells, we used 5 × 5 × 5 for the UO₂ unit cell; 2 × 5 × 5 for the 48-atom supercells (fewer divisions along the elongated direction and more along the shorter ones); and 2 × 2 × 2 for the 96-, 144-, and 216-atom supercells. A K-point convergence test on the 2 × 2 × 2 supercell (SI Fig. S1) shows that a 2 × 2 × 2 mesh meets our accuracy targets. With the increase of the supercell size, the Brillouin zone shrinks, so this mesh provides an equal or higher k-point sampling density for the larger cells. For all self-consistent field (SCF) calculations, the energy convergence threshold was set to 2 × 10⁻⁶ eV.

NEP configuration

NEP encodes the relative positions and interactions of each atom with its neighbors into compact descriptors and maps these local environments to per-atom energies via a single-hidden-layer feed-forward neural network.³⁵ This compact architecture enables efficient regression while retaining sufficient flexibility. The total potential energy is obtained by summing the energy U_i of each particle i, which is given by:where N_neu is the number of hidden neurons and N_des is the dimension of the input features. ω⁽⁰⁾ and ω⁽¹⁾ are the weights from the input to the hidden layer and from the hidden to the output layer, respectively. The hyperbolic tangent is used as the activation function to capture complex non-linear relationships. The input descriptor vector qⁱ is invariant to global translation, rotation, and permutation of identical atoms. b⁽⁰⁾ and b⁽¹⁾ are the biases of the hidden and output layers, respectively.

A clear definition of the loss function is critical during the optimization of model parameters. In NEP training, the loss is primarily based on the root mean square error (RMSE) between model predictions and reference values from first-principles calculations. To enhance interpolation and extrapolation capabilities, appropriate regularization is included. The loss function is defined as:

where λ_e, λ_f, λ_v, λ₁, and λ₂ are the weights for energy, atomic forces, virial tensor, L1 regularization, and L2 regularization, respectively. N_par, N_str, and N denote the numbers of trainable parameters, training structures, and atoms, respectively. z is an abstract vector composed of the trainable parameters. U^NEP(n,z) and W^NEP_μν(n,z) are the energy and virial tensor of structure n, and F^NEP_i(z) is the force vector on atom i. The corresponding reference quantities are U^ref(n), F^ref_i, and W^ref_μν(n). A more comprehensive description of NEP can be found in ref. 30.

Based on the preliminary tests of key hyperparameters (see SI Table S2), we finalized the NEP training configuration. The cutoff radii for the radial and angular descriptors were set to 6 Å and 5 Å, respectively. The neural network hidden layer contained 50 neurons. During training, the loss weights on energy, atomic forces, and virial tensor terms were set to 1, 1, and 0.1, respectively. Both L1 and L2 regularization were applied to enhance generalization. Other training hyperparameters are listed in SI Table S3. The overall NEP training workflow is shown in Fig. 1.

Fig. 1 Construction of the initial UO₂ dataset, the active learning workflow of NEP model training, and the structural types along with their counts in the final dataset.

To assess dataset representativeness, we analyzed the final dataset using scikit-learn.⁴⁴ As shown in Fig. 2a, the density distribution is approximately the Gaussian profile, indicating good coverage of the structural space. Fig. 2b shows a principal component analysis (PCA) projection of the UO₂ structural descriptors. This analysis demonstrates that active learning effectively broadened structural sampling, enabling exploration of a wider space phase. Finally, after 200 000 steps of full-batch training, we obtained the NEP model applicable to bulk UO₂ across a temperature range of approximately 0–2500 K.

Fig. 2 Statistical analysis of the final training dataset. (a) Sample density distribution of the training dataset. The vertical axis indicates the number of configurations. (b) Projected schematic of structural descriptors. Each point represents an individual structure, with color encoding its energy. The region enclosed by the solid blue line denotes the artificially constructed initial configurations (including lattice distortions and atomic perturbations), marked with squares; the region enclosed by the solid red line denotes configurations obtained from MD simulations and active learning, marked with stars.

Molecular dynamics simulations

All molecular dynamics simulations were performed with the GPUMD package using a 0.5 fs time step at a target pressure of 0 GPa. Unless otherwise specified, we used an 8 × 8 × 8 supercell (6144 atoms; a = b = c = 42.98 Å) of the conventional fluorite UO₂ unit cell. NPT (isothermal–isobaric) simulations employed the stochastic cell rescaling (SCR)⁴⁵ barostat, and NVT (canonical) simulations employed the Bussi–Donadio–Parrinello (BDP)⁴⁶ thermostat.

For structural properties, radial distribution functions (RDFs) were obtained after equilibrating the system for 100 ps in NPT–SCR and sampling for an additional 100 ps in NVT–BDP; RDFs were computed over 0–7 Å. Thermomechanical behavior was characterized by elastic constants, evaluated by applying a small strain of 0.01 to configurations that had first relaxed for 200 ps in NPT–SCR and then equilibrated for 400 ps in the NVE (microcanonical) ensemble. Thermal responses were quantified from NPT–SCR thermal-expansion simulations carried out from 200 to 2100 K in 100 K increments. At each temperature, we ran 400 ps of MD (200 ps equilibration and 200 ps production) and recorded thermodynamic quantities every 1 ps to extract the heat capacity and the linear thermal expansion coefficient.

For thermal transport, the thermal conductivity was cross validated by EMD, HNEMD, and NEMD. We assessed finite-size effects and the dependence on the magnitude and direction of the driving force (SI Fig. S2–S4). In EMD, a 200 ps NPT–SCR equilibration preceded a 1 ns NVT–BDP production run; the heat flux was sampled every 20 fs with a 100 ps correlation time. In HNEMD, we performed a 2 ns NVT–BDP production run after 200 ps of NPT–SCR equilibration, applying a uniform per-atom driving force of 2 × 10⁻⁴ Å⁻¹ along the y-direction and recording running averages of the thermal conductivity every 2 ps. In NEMD, heat transport along y was modeled using 8 × 26 × 8 (19 968 atoms; a = c = 42.98 Å, b = 139.68 Å) and 8 × 30 × 8 (23 040 atoms; a = c = 42.98 Å, b = 161.17 Å) supercells; the system was partitioned into eight groups along y (SI Fig. S5). After 1 ns of NVT–BDP equilibration, Langevin heat baths⁴⁷ were applied to the first and last groups as a heat source and cold sink, and the system was propagated for 2 ns. For both HNEMD and NEMD, spectral decomposition along y was performed with a correlation-time window of ±2 ps.

Results and discussion

Model training

Using the fourth-generation NEP framework as implemented in GPUMD, we performed full-batch training for 200 000 steps on the training dataset comprising 1918 structures. As shown in Fig. 3a, the loss functions for energies, atomic forces, and virial tensors converged during training, while the L1 and L2 regularization terms stabilized, indicating robust convergence of the NEP model. As presented in Fig. 3b–d, NEP predictions are in close agreement with DFT reference values, with RMSEs of 3.16 meV per atom for energies, 115.02 meV Å⁻¹ for forces, and 22.16 meV per atom for virial tensors. Across these metrics, NEP outperforms previously reported HDNNP,²⁰ SNAP²⁰ and NNP²² potentials. In addition, we assembled an independent test dataset of 200 structures for comparative evaluation (see SI Fig. S6). The RMSEs on this test dataset are likewise low, with no signs of overfitting, preliminarily indicating that the model generalizes well to new samples drawn from the same distribution as the training dataset. Overall, the NEP model trained in this study exhibits high predictive accuracy and robustness.

Fig. 3 Training results for the NEP model. (a) Convergence of the loss functions for energies, forces, and virial tensors, and of the L1 and L2 regularization terms during training. (b)–(d) Comparisons between NEP predictions and DFT results for (b) energies, (c) forces, and (d) virial tensors. The values in the subplots are root mean square errors.

Model validation

To systematically evaluate the reliability of the developed UO₂–NEP model, we selected three metrics that are highly sensitive to the accuracy of potential and have clear physical significance, and used them to benchmark the NEP model against DFT. These metrics include the zero-temperature equation of state (EOS), single-crystal elastic constants, and phonon dispersion relations across the full Brillouin zone. The EOS describes the response of materials under variations in pressure, volume, and temperature, provides a fundamental basis for determining the equilibrium lattice constant, quantifying compressibility, analyzing phase transition behavior, and assessing other thermodynamic properties. We sampled 12 volume points within the lattice scaling factor range of 0.94 to 1.05 and performed static total energy calculations on the corresponding structures to obtain EOS. Fig. 4 shows that the energy–volume (E–V) and pressure–volume (P–V) curves from NEP and DFT exhibit excellent agreement near the equilibrium volume. This indicates that the NEP model accurately captures the E–V relationship and pressure response of UO₂, providing a reliable foundation for subsequent simulations of thermodynamic properties.

Fig. 4 Energy–volume and pressure–volume curves for UO₂ at 0 K computed with NEP (blue) and DFT (yellow).

Elastic constants are critical indicators for evaluating the mechanical fidelity of an interatomic potential, as they are directly related to the mechanical performance and structural stability of materials. In this study, the elastic constants were calculated from the linear stress–strain constitutive relation. Within the small stress–strain regime, stress and strain obey Hooke's law:

where σ_i and ε_j denote the stress and strain components, and c_ij are the components of the elastic stiffness tensor in Voigt notation. To compute the elastic constants, a series of small strains ranging from −0.01 to 0.01 in increments of 0.005 were applied to the relaxed UO₂ unit cell, and the corresponding stresses were calculated. For cubic UO₂, the bulk modulus B, shear modulus G, Young's modulus E, and Poisson's ratio ν were estimated using the Voigt–Reuss–Hill (VRH)⁴⁸ averaging scheme. The Hill averages are reported as the final values (see Table 1). The results show that NEP predicts the UO₂ elastic constants C11, C12, and C44 with deviations relative to DFT of −1.16%, 7.43%, and 5.32%, respectively, while the relative errors in the remaining moduli (B, G, E, ν) are all within 2%. The CRG¹⁹ potential reproduces the experimental data well, with all elastic constants and moduli within 3% of experimental values. By contrast, both the NEP and MLIP-DFT²¹ models generally overestimate shear-related quantities relative to experimental data. Compared with the MLIP-DFT²¹ model, the NEP shows a more balanced performance in terms of overall accuracy and consistency across different moduli. Although its agreement with DFT for C44 and the bulk modulus is slightly inferior to that of MLIP-DFT,²¹ it performs better for G, E, and ν, and its bulk modulus is in closer agreement with the experiment value. These results suggest that the NEP model reliably captures the elastic response of UO₂, with particularly high fidelity for bulk compressibility.

Table 1 Comparison of elastic constants (C11, C12, C44), bulk modulus (B), shear modulus (G), Young's modulus (E), and Poisson's ratio (ν) for UO₂ at 0 K from NEP, DFT, the MLIP-DFT²¹ model, and the CRG¹⁹ empirical potential, alongside experimental⁴⁹ values at 25 °C

	NEP	DFT	MLIP-DFT	CRG	Exp.
C11 (GPa)	400.5	405.2	389.5	406.3	395.0
C12 (GPa)	118.6	110.4	121.2	124.7	121.0
C44 (GPa)	81.1	77.0	78.0	63.9	64.1
B (GPa)	212.5	208.6	207.8	218.6	212.3
G (GPa)	101.4	100.1	97.1	88.1	87.4
E (GPa)	262.4	259.0	252.0	233.0	230.6
ν	0.294	0.293	0.298	0.322	0.319

---

Phonon dispersion relations are a key tool for investigating lattice dynamical behavior and are crucial for understanding thermal transport. To evaluate the predictive performance of the NEP model, we computed the phonon dispersion relations of UO₂ using both the NEP model and DFT. After high-precision optimization of the UO₂ unit cell, the optimized structure was expanded to a 2 × 2 × 2 supercell for both calculations. NEP-based phonon dispersions were obtained in GPUMD with the finite-displacement method using a displacement of 0.01 Å. For the DFT reference, it was processed with PHONOPY⁵⁰ using the same displacement to construct the force-constant matrix. As shown in Fig. 5, neither approach yields imaginary frequencies throughout the Brillouin zone, indicating that the UO₂ structure is dynamically stable. In the low-frequency region, the acoustic branches show excellent agreement with the DFT results, demonstrating that the NEP model accurately captures the behavior of acoustic phonons, which dominate lattice thermal transport.²

Fig. 5 Phonon dispersion relations for UO₂ at 0 K obtained from the NEP model (solid blue lines) and DFT (yellow circles).

Radial distribution function

The radial distribution function (RDF) describes the probability density of the relative distances between any two particles in a many-particle system, thereby providing a direct measure of local structural order and the coordination environment. For a system of volume V containing N_a atoms of type a and N_b atoms of type b, the RDF g_ab(r) for the a–b pair is defined as:where dN_ab(r,dr) is the number of a–b pairs with separations in the interval [r, r + dr]. Fig. 6 shows that the peak heights of all pairs decrease with increasing temperature, primarily due to enhanced thermal motion leading to greater structural disorder. The positions of the first peak for U–U and O–O pairs shift slightly to larger r as temperature rises, indicating an increase in the average interatomic distance, which reflects thermal expansion. In contrast, the first peak of the U–O pair shifts slightly to smaller r at higher temperatures, suggesting a minor contraction of the average bond length and a concomitant relative strength of the U–O bond. These trends are consistent with the NNP²² model.

Fig. 6 Partial (U–U, O–O, and U–O) and total radial distribution functions (RDFs) over the temperature range of 300–2500 K.

Elastic constants

We systematically calculated the three independent isothermal elastic constants, C11, C12, and C44, over the range of 300–2000 K (see Fig. 7). Overall, the elastic constants predicted by different force-field models show deviations of varying magnitude from experimental values. These discrepancies primarily arise from: (I) inherent limitations of the interatomic potentials, including approximations to the potential energy surface and limited coverage of the training datasets; (II) intrinsic constraints of MD simulations, such as finite-size effects, implementation differences between the stress–strain and fluctuation approaches, and comparability between the computed isothermal elastic constants and the adiabatic elastic constants commonly measured experimentally; (III) technical challenges in experimental measurements, including uncertainty at high temperature, non-stoichiometry in UO_2±x, defects, and effects of polycrystalline; and (IV) the intrinsic complexity of the mechanical behavior with temperature and microstructural evolution.

Fig. 7 Comparison of the elastic constants C11, C12, and C44 as a function of temperature (300–2000 K). NEP (blue circles), NNP²² (yellow squares), SNAP²⁰ (green upward triangles), and HDNNP²⁰ (red downward triangles) are machine learning potentials; CRG¹⁹ (purple diamonds) is an empirical potential; experimental values from Hutchings⁵¹et al. are shown as brown crosses.

Specifically, for C11, the HDNNP²⁰ and SNAP²⁰ models systematically underestimate the experimental values over most of the temperature range, whereas NNP,²² CRG,¹⁹ and the NEP models developed in this study tend to overestimate slightly. Notably, the experimentally reported scatter in C11 at 300 K is about 60 GPa, indicating substantial uncertainty in the experimental reference data, thereby complicating model assessment. For C12, except for a slight overestimation by the NNP²² model at 300 K, the other models generally lie below the experimental values. All models predict a monotonic decrease in C12 with increasing temperature, in marked contrast to the nearly constant experimental trend of around 120 GPa. For C44, the CRG,¹⁹ HDNNP,²⁰ and SNAP²⁰ models yield values lower than experimental values, whereas the NNP²² and NEP models show closer agreement with the experimental data. Overall, the NEP model exhibits smaller deviations from experimental data for both C11 and C44, showing a potential advantage in characterizing the elastic properties of UO₂.

Heat capacity and the linear thermal expansion coefficient

For three-dimensional isotropic systems, the linear thermal expansion coefficient α_L at constant pressure is approximated here by a central finite difference:

Similarly, the isobaric molar heat capacity C_p is given by:

where T_k (2 ≤ k ≤ 19) denotes the target temperature of segment k, ΔT is the interval between adjacent temperature segments, n is the number of moles in the system, L(T_k) is the average length along any crystal axis, and H(T_k) is the enthalpy of the system. Mean values and standard errors (SE) of these physical quantities were computed from the final 200 production data points in each simulation segment. As shown in Fig. 8a, the C_P values predicted by NEP are consistent with those obtained from the SNAP²⁰ model. Above 1500 K, both the NEP and the MLIP-DFT²¹ models exhibit a pronounced upward trend. Over the range of 300–2000 K, the NEP results are generally consistent with experimental measurements. Fig. 8b indicates that the linear thermal expansion coefficients computed using NEP are in close agreement with those derived from the Basak¹⁶ potential up to 1000 K, beyond which the NEP predicted coefficients increase rapidly with temperature.

Fig. 8 Temperature dependence of (a) the isobaric molar heat capacity and (b) the linear thermal expansion coefficient. The black circles with error bars represent the NEP results from this work. Crosses denote experimental reference values, taken from Ronchi⁵ and Martin et al.⁵² Data for other potentials are taken from ref. 16,19–22.

Thermal conductivity

Thermal conductivity is a critical parameter that characterizes the heat transfer capability and is essential for the design and optimization of nuclear fuels. In this work, we employed EMD, HNEMD, and NEMD to compute the thermal conductivity of UO₂ over a range of temperatures.

In the EMD approach, the thermal conductivity is obtained from the Green–Kubo relation via the heat current autocorrelation function (HCACF) of an equilibrium system. The Green–Kubo expression reads:

where κ_μν(t)(μ,ν = x, y, z) is the thermal conductivity tensor, k_B is the Boltzmann constant, T and V are the temperature and volume, and 〈J_μ(0)J_ν(t′)〉 denotes the HCACF. When t is sufficiently long, the running integral κ_μν(t) converges to the thermal conductivity tensor of the system. In GPUMD, an efficient optimized algorithm³¹ is implemented for heat flux calculations, which derives the heat current from per-atom virial tensors, and GPU acceleration is used to evaluate the HCACF and its time integral, thereby avoiding redundant computations. The HNEMD method is suitable for studying heat transport in homogeneous systems subject to an external driving field. By applying an additional driving force, a nonequilibrium steady state is established, yielding the heat flux J:

The thermal conductivity is then given by:

where F_e is the homogeneous driving force. It must be small enough to keep the system within the linear-response regime and large enough to retain a sufficiently large signal-to-noise ratio. Within the sufficiently long simulation times and suitable F_e, this expression is equivalent to the Green–Kubo results. Because the heat current is the time derivative of the energy moment, a finite F_e continuously injects power into the system. Consequently, HNEMD must be coupled to an unbiased thermostat to remove this additional power and maintain a steady state.³² In NEMD, a steady-state temperature gradient is established by imposing heat source and cold sink regions at opposite ends of the system, thereby driving heat flow. This approach is formally consistent with Fourier's law, and the thermal conductivity at thermal transport channel L is expressed as:where G(L) is the thermal conductance per unit area, Q is the steady-state heat current, S is the cross-sectional area perpendicular to the transport direction, and ΔT/L denotes the end-to-end temperature gradient. Since the end-to-end temperature difference ΔT includes contributions from the non-linear temperature regions near the heat source and cold sink, as well as the temperature drop due to the contact thermal resistance (Kapitza jump), the quantity κ(L) defined in this way is not an intrinsic material property but an apparent thermal conductivity that depends on the finite-size L. Further theoretical and implementation details can be found in ref. 33.

To efficiently obtain the apparent thermal conductivity κ(L) at arbitrary lengths, we take the Fourier transforms of the virial-velocity cross-correlation function K(t) in both HNEMD and NEMD, thereby establishing a unified spectral decomposition framework. Within this framework, the spectral thermal conductivity κ(ω) in HNEMD is defined as:

Correspondingly, the spectral thermal conductance G(ω) along the same transport direction in NEMD is given by:

The spectral mean free path λ(ω) is then defined as:

The length-dependent thermal conductivity κ(L) is then obtained from the above spectral quantities:

In this study, we computed the thermal conductivity at 500 K using EMD, HNEMD, and NEMD. The results are shown in Fig. 9. For the EMD results shown in Fig. 9a, because the method relies on integrating the HCACF, the long tail exhibits a low signal-to-noise ratio and is prone to drift, necessitating multiple independent simulations to reduce uncertainty. We conducted 30 independent runs reducing the uncertainty to about 7%. In contrast, the HNEMD in Fig. 9b directly drives a steady-state heat flux via a perturbative field, yielding stronger signals and faster convergence. Only three independent simulations were required to obtain stable estimates, with markedly lower statistical uncertainty than in EMD. NEMD requires large supercells to establish a stable temperature gradient, leading to the highest computational cost per simulation. SI Fig. S7 verifies satisfactory group thermostatting; the spectral thermal conductivity from HNEMD and the spectral thermal conductance from NEMD are shown in SI Fig. S8 and S9, respectively. Fig. 9c shows that the spectral phonon mean free path in UO₂ at 500 K spans roughly one order of magnitude, from approximately 50 nm at high frequencies to approximately 500 nm in the low-frequency limit. Correspondingly, in Fig. 9d, κ(L) increases with L and approaches a constant for L exceeding about 5 µm. This asymptotic value can be regarded as bulk (macroscopic) thermal conductivity. Consistent with this, Watanabe et al. reported a pronounced finite-size dependence of thermal conductivity on characteristic length.⁵³ This behavior originates from the temperature-dependent distribution of the phonon mean free path: at high temperatures, enhanced Umklapp scattering shortens the mean free path, making κ less sensitive to the geometric size, whereas at low temperatures the fraction of long-mean-free-path phonons increases, allowing them to propagate over longer distances with fewer scattering events.

Fig. 9 (a) and (b) Mean values and uncertainties in the thermal conductivity of UO₂ at 500 K obtained from EMD and HNEMD, respectively; (c) spectral phonon mean free path; and (d) apparent thermal conductivity derived from HNEMD, compared with the NEMD results.

Given the comparative efficiency of the three methods assessed above, we adopted the HNEMD approach, which possesses faster convergence and lower uncertainty, to simulate the thermal conductivity of UO₂ over 300–1500 K. Ronchi et al.⁵ noted that the dominant thermal transport mechanisms vary with temperature: lattice contributions prevail at low temperatures, whereas non-lattice channels such as electronic contributions become markedly enhanced above 1300 K. Since classical MD explicitly describes phonons and diffusive carriers (e.g., ions, defects) but excludes electronic contributions, we restricted our simulations to temperatures below 1500 K. As shown in Fig. 10, the NEP model is in overall agreement with the experimental results of Ronchi et al.⁵ and with the EIP predictions of Zhou et al.⁵⁴ above 800 K. Below 800 K, the NEP slightly overestimates the experimental values yet remains lower than the results obtained from the CRG¹⁹ empirical potential and the modified Busker potential of Liu et al.⁵⁵ (fitted via the Callaway model with spin phonon relaxation corrections), and it aligns closely with the PBE+U calculations reported by Zhou et al.⁵⁴ These discrepancies mainly stem from the use of an ideal, defect-free model, which omits impurities, defects, and electronic contributions present in experimental specimens. We note that electronic thermal conduction is non-negligible in experiments within the 1300–1500 K range, which partly explains why our predictions are slightly lower than the measurements in this range. Overall, the NEP model combined with the HNEMD reproduces experimental trends at intermediate temperatures and yields predictions consistent with first-principles calculations at lower temperatures, demonstrating its utility for characterizing the thermal transport properties of nuclear fuels.

Fig. 10 Thermal conductivity of UO₂ over 300–1500 K. The solid lines represent the results obtained from the CRG¹⁹ empirical potential (blue), PBE+U⁵⁴ calculations reported by Zhou et al. (yellow), the EIP⁵⁴ of Zhou et al. (green), and the modified Busker⁵⁵ potential of Liu et al. (red). Experimental data points are from Fink et al.⁶ (blue crosses), Ronchi et al.⁵ (yellow squares), and Lucuta et al.⁴ (green plus signs). NEP predictions are shown as purple circles, with error bars denoting uncertainties.

Conclusions

The actinide oxide UO₂ is the primary fuel in commercial nuclear reactors, and its thermal transport properties are pivotal to reactor-core thermal management and safety margins. In this work, we developed an NEP machine learning potential approaching DFT-level accuracy and combined with the energy-conserving Hardy-type heat flux framework in GPUMD, and we systematically evaluated and cross-validated the thermal conductivity of UO₂.

Initially, by manually generating an initial dataset and adopting an active learning strategy of “training–sampling–labeling”, we constructed a training dataset that covers both static and dynamic configurations with a balanced density distribution. After iterative refinements, we obtained an NEP applicable for UO₂ over 0–2500 K, which agrees closely with DFT for predicted energies, forces, and virials. Further validation confirmed that the equation of state derived from NEP closely matches that of DFT near the equilibrium volume. The predicted elastic constants C11, C12, and C44 deviate from DFT reference values by −1.16%, 7.43%, and 5.32%, respectively. Phonon dispersion relations show no imaginary modes across the Brillouin zone, and the acoustic branches align well with DFT. Overall, the NEP exhibits high accuracy, strong robustness, and reliable capability for characterizing the fundamental physical properties of UO₂.

Subsequently, we further examined the temperature dependence of various thermophysical properties. The radial distribution functions indicate that the first-neighbor distances for U–U and O–O pairs increase with increasing temperature, while the average U–O bond length decreases slightly. For elastic constants, the NEP predictions for C11 and C44 are closer to experimental data than those from other MLIPs, and C12 decreases monotonically with temperature. The heat capacity shows good agreement with experimental data in the range of 300–2000 K. The linear thermal expansion coefficient agrees with the Basak potential at T ≤ 1000 K and increases rapidly at higher temperatures. These results underscore the strong applicability and promising potential of the NEP for characterizing the thermophysical properties of UO₂.

Finally, we systematically evaluated the EMD, NEMD, and HNEMD methods for thermal conductivity calculations and cross-validated the consistency of the NEP predictions across methods. Among them, the HNEMD delivers the highest computational efficiency, characterized by a high signal-to-noise ratio, low uncertainty, and rapid convergence, outperforming EMD (which requires multiple independent repetitions) and NEMD (with substantial length-scale requirements). At 500 K, the spectral phonon mean free path distribution spans roughly one order of magnitude, and the thermal conductivity reaches the bulk limit once the transport channel length exceeds about 5 µm. Based on HNEMD, the NEP-predicted thermal conductivity from 300 to 1500 K shows good agreement with experimental data above 800 K and aligns closely with the PBE+U results from Zhou et al. below 800 K. Thus, the NEP combined with HNEMD reliably captures the intermediate-temperature thermal transport behavior of UO₂ while balancing efficiency and accuracy and yields predictions near first-principles precision in the low-temperature regime.

This work demonstrates the applicability and utility of the NEP for modeling thermal transport in nuclear fuels and provides methodological support for multiscale integration. Beyond UO₂, the same framework is applicable to other actinide oxide fuels, as well as to non-stoichiometric compositions, interfaces and grain boundaries, and irradiation-damaged microstructures that govern heat conduction. Such applications can help assess transferability across chemistries and operating conditions and broaden the use of NEP-based simulations in fuel-performance modeling.

Author contributions

Junying Zhong: writing – original draft, writing – review and editing, investigation, methodology, software, visualization, data curation, formal analysis, and validation. Lei Zhang: investigation, conceptualization, resources, and supervision. Tao Bo: writing – review and editing, data curation, methodology, investigation, conceptualization, software, resources, formal analysis, validation, supervision, project administration, and funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this study are included as part of the Supplementary information (SI). The SI also provides the GPUMD input files for all MD simulations (RAR archive with nep.txt, model.xyz, a README, and per-property folders each containing run.in). See DOI: https://doi.org/10.1039/d5cp03760c.

Acknowledgements

We acknowledge financial support from the National Natural Science Foundation of China (No. 12575387, U2167218 and U2267222) and the Yongjiang Talent Introduction Programme (2021A-161-G).

Notes and references

1. X. Yang, J. Tiwari and T. Feng, Mater. Today Phys., 2022, 24, 100689.
2. D. H. Hurley, A. El-Azab, M. S. Bryan, M. W. D. Cooper, C. A. Dennett, K. Gofryk, L. He, M. Khafizov, G. H. Lander and M. E. Manley, et al. , Chem. Rev., 2022, 122, 3711–3762.
3. J. H. Harding and D. G. Martin, J. Nucl. Mater., 1989, 166, 223–226.
4. P. G. Lucuta, H. Matzke and I. J. Hastings, J. Nucl. Mater., 1996, 232, 166–180.
5. C. Ronchi, M. Sheindlin, M. Musella and G. J. Hyland, J. Appl. Phys., 1999, 85, 776–789.
6. J. K. Fink, J. Nucl. Mater., 2000, 279, 1–18.
7. J. J. Carbajo, G. L. Yoder, S. G. Popov and V. K. Ivanov, J. Nucl. Mater., 2001, 299, 181–198.
8. M. Sanati, R. C. Albers, T. Lookman and A. Saxena, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 014116.
9. Y. Yun, D. Legut and P. M. Oppeneer, J. Nucl. Mater., 2012, 426, 109–114.
10. B.-T. Wang, P. Zhang, R. Lizárraga, I. Di Marco and O. Eriksson, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88, 104107.
11. B.-T. Wang, J.-J. Zheng, X. Qu, W.-D. Li and P. Zhang, J. Alloys Compd., 2015, 628, 267–271.
12. T. P. Kaloni, N. Onder, J. Pencer and E. Torres, Ann. Nucl. Energy, 2020, 144, 107511.
13. S. Sheykhi and M. Payami, Iran. J. Sci. Technol., Trans. A: Sci., 2020, 44, 1585–1593.
14. E. Torres, I. Cheik Njifon, T. P. Kaloni and J. Pencer, Comput. Mater. Sci., 2020, 177, 109594.
15. S. T. Murphy, M. J. D. Rushton and R. W. Grimes, Prog. Nucl. Energy, 2014, 72, 27–32.
16. C. B. Basak, A. K. Sengupta and H. S. Kamath, J. Alloys Compd., 2003, 360, 210–216.
17. N.-D. Morelon, D. Ghaleb, J.-M. Delaye and L. Van Brutzel, Philos. Mag., 2003, 83, 1533–1555.
18. E. Yakub, C. Ronchi and D. Staicu, J. Chem. Phys., 2007, 127, 094508.
19. M. W. D. Cooper, M. J. D. Rushton and R. W. Grimes, J. Phys.: Condens. Matter, 2014, 26, 105401.
20. E. T. Dubois, J. Tranchida, J. Bouchet and J.-B. Maillet, Phys. Rev. Mater., 2024, 8, 025402.
21. E. Stippell, L. Alzate-Vargas, K. N. Subedi, R. M. Tutchton, M. W. D. Cooper, S. Tretiak, T. Gibson and R. A. Messerly, Artif. Intell. Chem., 2024, 2, 100042.
22. K. Konashi, N. Kato, K. Mori and K. Kurosaki, J. Nucl. Mater., 2025, 607, 155660.
23. K. Govers, S. Lemehov, M. Hou and M. Verwerft, J. Nucl. Mater., 2008, 376, 66–77.
24. S. Nichenko and D. Staicu, J. Nucl. Mater., 2014, 454, 315–322.
25. C. I. Maxwell and J. Pencer, Ann. Nucl. Energy, 2019, 131, 317–324.
26. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in’t Veld, A. Kohlmeyer, S. G. Moore and T. D. Nguyen, et al. , Comput. Phys. Commun., 2022, 271, 108171.
27. A. P. Thompson, S. J. Plimpton and W. Mattson, J. Chem. Phys., 2009, 131, 154107.
28. P. Boone, H. Babaei and C. E. Wilmer, J. Chem. Theory Comput., 2019, 15, 5579–5587.
29. S. T. Tai, C. Wang, R. Cheng and Y. Chen, J. Chem. Theory Comput., 2025, 21, 3649–3657.
30. Z. Fan, Y. Wang, P. Ying, K. Song, J. Wang, Y. Wang, Z. Zeng, K. Xu, E. Lindgren and J. M. Rahm, et al. , J. Chem. Phys., 2022, 157, 114801.
31. Z. Fan, L. F. C. Pereira, H.-Q. Wang, J.-C. Zheng, D. Donadio and A. Harju, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 094301.
32. Z. Fan, H. Dong, A. Harju and T. Ala-Nissila, Phys. Rev. B, 2019, 99, 064308.
33. Z. Li, S. Xiong, C. Sievers, Y. Hu, Z. Fan, N. Wei, H. Bao, S. Chen, D. Donadio and T. Ala-Nissila, J. Chem. Phys., 2019, 151, 234105.
34. Z. Fan, Z. Zeng, C. Zhang, Y. Wang, K. Song, H. Dong, Y. Chen and T. Ala-Nissila, Phys. Rev. B, 2021, 104, 104309.
35. Z. Fan, J. Phys.: Condens. Matter, 2022, 34, 125902.
36. K. Song, R. Zhao, J. Liu, Y. Wang, E. Lindgren, Y. Wang, S. Chen, K. Xu, T. Liang and P. Ying, et al. , Nat. Commun., 2024, 15, 10208.
37. C. Chen, Y. Li, R. Zhao, Z. Liu, Z. Fan, G. Tang and Z. Wang, Comput. Phys. Commun., 2025, 317, 109859.
38. T. Liang, K. Xu, E. Lindgren, Z. Chen, R. Zhao, J. Liu, E. Berger, B. Tang, B. Zhang, Y. Wang and et al., arXiv, 2025, preprint, arxiv:2504.21286 DOI:10.48550/arXiv.2504.21286.
39. J. Hafner, J. Comput. Chem., 2008, 29, 2044–2078.
40. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
41. A. Dal Corso, A. Pasquarello, A. Baldereschi and R. Car, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 53, 1180–1185.
42. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
43. H. J. Monkhorst and J. D. Pack, Phys. Rev. B, 1976, 13, 5188–5192.
44. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss and V. Dubourg, et al. , J. Mach. Learn. Res., 2011, 12, 2825–2830.
45. M. Bernetti and G. Bussi, J. Chem. Phys., 2020, 153, 114107.
46. G. Bussi, D. Donadio and M. Parrinello, J. Chem. Phys., 2007, 126, 014101.
47. G. Bussi and M. Parrinello, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 056707.
48. R. Hill, Proc. Phys. Soc., London, Sect. A, 1952, 65, 349–354.
49. J. B. Wachtman Jr., M. L. Wheat, H. J. Anderson and J. L. Bates, J. Nucl. Mater., 1965, 16, 39–41.
50. A. Togo, J. Phys. Soc. Jpn., 2023, 92, 012001.
51. M. T. Hutchings, J. Chem. Soc., Faraday Trans. 2, 1987, 83, 1083–1103.
52. D. G. Martin, J. Nucl. Mater., 1988, 152, 94–101.
53. T. Watanabe, S. B. Sinnott, J. S. Tulenko, R. W. Grimes, P. K. Schelling and S. R. Phillpot, J. Nucl. Mater., 2008, 375, 388–396.
54. S. Zhou, C. Jiang, E. Xiao, S. Bandi, M. W. D. Cooper, M. Jin, D. H. Hurley, M. Khafizov and C. A. Marianetti, J. Phys.: Condens. Matter, 2025, 37, 255901.
55. X.-Y. Liu, M. W. D. Cooper, K. J. McClellan, J. C. Lashley, D. D. Byler, B. D. C. Bell, R. W. Grimes, C. R. Stanek and D. A. Andersson, Phys. Rev. Appl., 2016, 6, 044015.

---