Warm dense matter (WDM) is a state of matter lying between condensed matter and plasma and consisting of strongly coupled ions and partially degenerate electrons. WDM exists in the interiors of giant planets^1,2 and in the crusts of white dwarf and neutron stars,^3,4 and it can be generated through laboratory experiments with diamond anvil cells,⁵ gas guns,^6,7 and high-power lasers.^7,8 WDM also plays an important role in inertial confinement fusion (ICF).⁹ In these contexts, it is crucial to understand the properties of WDM such as its equation of state and its optical and transport properties. Owing to the lack of experimental data, quantum-mechanics-based simulation methods such as Kohn–Sham density functional theory (KSDFT),^10,11 orbital-free density functional theory (OFDFT),¹² and path-integral Monte Carlo^13–16 have emerged as ideal tools to study WDM.^16–20 Thermal conductivity is one of the most important properties of WDM. For example, it has been widely studied in the context of models of planetary interiors.^21–23 It is a key parameter governing the growth of hydrodynamic instabilities that arise with various capsule designs for the National Ignition Facility.²⁴ It also plays a key role in simulations of the interactions between lasers and metal targets.²⁵ The thermal conductivity has a significant influence on predictions of ICF implosions in hydrodynamic simulations.²⁶

The thermal conductivity includes both electronic and ionic contributions. The Kubo–Greenwood (KG) formula^27,28 has been widely used in studies of the electronic thermal conductivity κ_e of liquid metals and WDM.^{18,26,29–35} Typically, the κ_e value is an average of results from several atomic configurations, which are selected from first-principles molecular dynamics (FPMD) simulations. However, it is computationally expensive to perform FPMD simulations for large systems, especially for WDM when the temperatures are high.^19,33,36,37 Besides, it has been shown that size effects may substantially affect κ_e.^30,33,38 For computations of the ionic thermal conductivity κ_I, we focus on the Green–Kubo (GK) formula,^27,39,40 where κ_I is expressed in terms of the heat flux autocorrelation function, which requires knowledge of the energy and the virial tensor of each atom from the molecular dynamics trajectory. Typically, this formula is used with empirical force fields,^41–43 since traditional DFT methods cannot yield an explicit energy for each atom; note that ab initio GK formulas for evaluation of κ_I have been proposed in some recent work.^44–47 Besides, a long molecular dynamics trajectory may be needed to evaluate κ_I, which poses another challenge for DFT simulations.

While current experimental techniques only measure the total thermal conductivity, first-principles methods have become ideal tools to yield electronic and ionic contributions separately. However, only a few studies adopting first-principles methods have examined transport properties by considering both contributions.^32,48,49 For some metals, the ionic contribution is not important,⁴⁸ but this is not the case for all metals, tungsten being one example.⁴⁹ In addition, first-principles methods are computationally expensive, especially when a large system or a long trajectory is considered. In this regard, there is a need for an efficient and accurate method that can yield both electronic and ionic thermal conductivities of materials. Here, using the KG and GK formulas, we propose a combined DFT and deep potential (DP) method for this purpose, taking warm dense Al as an example. The recently developed DP molecular dynamics (DPMD)^50–52 method learns first-principles data via deep neural networks and yields a highly accurate many-body potential to describe the interactions among atoms. Compared with traditional DFT, DPMD has a much higher efficiency while retaining the ab initio accuracy. In addition, linear-scaling DPMD can be parallelized to simulate hundreds of millions of atoms.^53,54 The DP method has been applied in a variety of contexts, such as crystallization of silicon,⁵⁵ high-entropy materials,⁵⁶ and isotope effects in liquid water,^57,58 as well as WDM.^59,60

In this work, we demonstrate that the DP method in conjunction with DFT can be used to obtain both electronic and ionic thermal conductivities of warm dense Al. We use KSDFT, OFDFT, and DPMD to study the electronic and ionic thermal conductivities of warm dense Al at temperatures of 0.5 eV, 1.0 eV, and 5.0 eV with a density of 2.7 g/cm³. The electronic thermal conductivity can be accurately computed via the KG method based on the DPMD trajectories instead of the FPMD trajectories. Importantly, we systematically investigate convergence issues with regard to the number of k-points, the number of atoms, the broadening parameter, the exchange-correlation functionals, and the pseudopotentials, together with their effects on the determination of the electronic thermal conductivity via DPMD simulations. Furthermore, the ionic thermal conductivity can also be obtained via DPMD simulations, and we study convergence for different sizes of systems and different lengths of trajectories.

The remainder of the paper is organized as follows. In Sec. II, we briefly introduce KSDFT, OFDFT, and DPMD, as well as the setups of the simulations. We also briefly introduce the KG and GK formulas that are used in this work. In Sec. III, we present the electronic thermal conductivities computed using each of the three methods. We then give a thorough discussion of convergence issues as they affect the electronic thermal conductivity. Finally, we present results for the ionic thermal conductivity of warm dense Al. We conclude our work in Sec. IV.

The ground-state total energy within the formalism of DFT^10,11 can be expressed as a functional dependence of the electron density. The Hohenberg–Kohn theorems¹⁰ point out that the electron density that minimizes the total energy is the ground-state density, while the energy is the ground-state energy.¹⁰ Depending on the way in which the electron kinetic energy is treated, there are two DFT methods, namely, KSDFT¹¹ and OFDFT.¹² The electron kinetic energy in KSDFT does not involve the electron density explicitly, but is evaluated from the ground-state wave functions of electrons obtained by self-consistent iterations. On the other hand, in OFDFT, the electron kinetic energy is given approximately as an explicit density functional,^61–65 which enables a direct search for the ground-state density and energy. As a result, KSDFT has higher accuracy, but OFDFT is several orders more efficient, especially for large systems. Mermin extended DFT to finite temperatures,⁶⁶ in an approach that has been widely used to describe electrons in WDM.^18,29

We run 64-atom Born–Oppenheimer molecular dynamics (BOMD) simulations with KSDFT using the Quantum ESPRESSO 5.0.2 package.⁶⁷ The Perdew–Burke–Ernzerhof (PBE) exchange-correlation (XC) functional⁶⁸ is used. The projector augmented-wave (PAW) potential^69,70 is used, with three valence electrons and a cutoff radius of 1.38 Å. The plane wave cutoff energy is set to 20 Ry for temperatures of 0.5 eV and 1.0 eV and to 30 Ry for 5.0 eV. We only use the gamma k-point. Periodic boundary conditions are adopted. The Andersen thermostat⁷¹ is used with the NVT ensemble, and the trajectory length is 10 ps, with a time step of 1.0 fs.

We also perform 108-atom BOMD simulations with OFDFT, using the PROFESS 3.0 package.⁷² Both PBE⁶⁸ and local density approximation (LDA)¹¹ XC functionals are used, together with the Wang–Teter (WT) kinetic energy density functional (KEDF)⁶⁵ and the Nosé–Hoover thermostat^73,74 in the NVT ensemble. We run OFDFT simulations for 10 ps, with a time step of 0.25 fs, and we set the energy cutoff to 900 eV. Periodic boundary conditions are adopted.

The DP method^50–52 learns the dependence of the total energy on the coordinates of the atoms in a system and constructs a deep neural network (DNN) model that predicts the potential energy and force of each atom. In the training process, the total potential energy E_tot is decomposed into the energies of the different atoms:Etot=∑iEi,(1)where E_i is the potential energy of atom i. In general, through the DNN, for each atom i, a mapping is established between E_i and the atomic coordinates of its neighboring atoms within a cutoff radius r_c. Specifically, the DNN model consists of an embedding network and a fitting network.⁵⁰ The embedding network imposes constraints on the atoms, such that the atomic coordinates obey symmetries under translation, rotation, and permutation. On the other hand, the fitting network maps the atomic coordinates from the embedding network to E_i. Next, the training data are prepared as the total energy and the forces acting on atoms are extracted from the FPMD trajectories. Finally, the parameters of the DNN model are optimized by minimizing a loss function. The resulting DNN-based model can be used to simulate a large system with an efficiency comparable to that when empirical force fields are used.

In this work, we adopt DNN-based models trained from either KSDFT or OFDFT trajectories with the DeePMD-kit package.⁵¹ With the purpose of studying size effects on electronic thermal conductivity, a series of cubic cells consisting of 16, 32, 64, 108, 216, and 256 atoms are simulated for 10 ps with a time step of 0.25 fs. In addition, we run systems of 16, 32, 64, 108, 256, 1024, 5488, 8192, 10 648, 16 384, 32 000, and 65 536 atoms to investigate the convergence of the ionic thermal conductivity. The length of these trajectories is 500 ps, with a time step of 0.25 fs. We employ the Nosé–Hoover thermostat^73,74 in the NVT ensemble. Periodic boundary conditions are adopted. All of the DPMD simulations are performed using the modified LAMMPS package.⁷⁵

The electronic thermal conductivity κ_e is calculated from the Onsager coefficients L_mn asκe=1e2TL22−L122L11,(2)where T is the temperature and e is the electron charge. The L_mn are obtained as the zero-frequency limits of the frequency-dependent Onsager coefficients L_mn(ω):Lmn=limω→0Lmn(ω).(3)We follow the approach to the KG formula adopted by Holst et al.,⁷⁶ according to which the L_mn(ω) are given byLmn(ω)=(−1)m+n2πe2ℏ23me2ωΩ×∑ijαkW(k)ϵik+ϵjk2−μm+n−2|〈Ψik|∇α|Ψjk〉|2×[f(ϵik)−f(ϵjk)]δ(ϵjk−ϵik−ℏω),(4)where m_e is the electron mass, Ω is the cell volume, W(k) is the weight of k-points in the Brillouin zone, μ is the chemical potential, Ψ_ik is the wave function of the ith band with eigenvalue ϵ_ik, and f is the Fermi–Dirac function. We use the chemical potential μ instead of the enthalpy per atom, which does not affect the results for the electronic thermal conductivity in one-component systems.⁷⁶ Note that here we use the momentum operator in Eq. (4) instead of the velocity operator, and this introduces additional approximations due to the use of nonlocal pseudopotentials.^{18,34,77–79} However, the nonlocal errors have been demonstrated to be about 3.3% for liquid Al,⁷⁸ and this is considered to be sufficiently small for the purposes of this work. In the future work, we will include the additional term caused by nonlocal pseudopotentials. We also define the frequency-dependent electronic thermal conductivity asκe(ω)=1e2TL22(ω)−L122(ω)L11(ω).(5)In practical applications of the KG method, the delta function in Eq. (4) needs to be broadened. We adopt a Gaussian function,²⁹ and the delta function then takes the formδ(E)=limΔE→012πΔEe−E2/2ΔE2.(6)Here, ΔE controls the full width at half maximum (FWHM, denoted by σ) of the Gaussian function through the relation σ ≈ 2.3548ΔE.

The KG method needs eigenvalues and wave functions computed from DFT solutions of given atomic configurations. In practice, we select 5–20 atomic configurations from the last 2 ps MD trajectories, with a time interval of 0.1 ps. We use both PBE and LDA XC functionals and the associated norm-conserving (NC) pseudopotentials to examine the influences of XC functionals and pseudopotentials on the resulting electronic thermal conductivity. We adopt two NC pseudopotentials for Al, which are referred to here as PP1 and PP2. The PP1 pseudopotential is generated with the optimized norm-conserving Vanderbilt pseudopotential method via the ONCVPSP package.^80,81 We use 11 valence electrons and a cutoff radius of 0.50 Å. Our calculations of κ_e in a 64-atom cell involve 1770, 2100, and 5625 bands at temperatures 0.5 eV, 1.0 eV, and 5.0 eV, respectively. The PP2 pseudopotential is generated through the PSlibrary package.⁸² We use the Troullier–Martins method,⁸³ and the cutoff radius is set to 1.38 Å. We choose three valence electrons for each atom. We select 720, 1100, and 4800 bands for calculations of κ_e at temperatures 0.5 eV, 1.0 eV, and 5.0 eV, respectively. The plane wave cutoff energies of both pseudopotentials are set to 50 Ry. Generally, the PP1 pseudopotential is used in most cases, with the PP2 pseudopotential only being used to compare the effects of different pseudopotentials on κ_e. Figure 1 shows the computed averaged κ_e(ω) with respect to different numbers of atomic configurations from 108-atom DPMD trajectories at 0.5 eV and 5.0 eV. The DP model is trained based on the 108-atom OFDFT trajectories using the PBE XC functional. We can see that 20 atomic configurations are enough to converge κ_e(ω). Additionally, κ_e(ω) converges more easily at the relatively low temperature of 0.5 eV. Therefore, we choose 20 atomic configurations for cells with 108 atoms or less at all the three temperatures. For systems with more than 108 atoms, we respectively select 5 and 10 atomic configurations for temperatures smaller than 5.0 eV and equal to 5.0 eV, unless otherwise specified.

In the DPMD method, the total potential energy of the system is decomposed onto each atom. In this regard, the ionic thermal conductivity can be calculated through the GK method⁴⁰ with the formulaκI=13ΩkBT2∫0+∞〈Jq(t)⋅Jq(0)〉dt,(7)where Ω is the cell volume, T is the temperature, k_B is the Boltzmann constant, 〈⋯〉 is the ensemble average, and t is time. J_q in Eq. (7) is the heat current of a one-component system in the center-of-mass frame and takes the formJq=∑i=1Nεivi−12∑i=1N∑j≠iN(vi⋅Fij)rij.(8)Here, v_i is the velocity of the ith atom, and ɛ_i is the energy of the ith atom including ionic kinetic energy and potential energy, where the potential energy is directly obtained by the DNN. F_ij is the force acting on the ith atom due to the presence of the jth atom, with r_ij defined as r_i − r_j. It is worth mentioning that Eq. (8) is valid for two-body potentials but yields an error of about 20% in J_q when a many-body potential (e.g., the deep potential) is adopted.⁸⁴ Importantly, as will be shown below, the ionic thermal conductivity of warm dense Al is at least two orders of magnitudes smaller than the electronic thermal conductivity. Therefore, we still consider Eq. (8) as a valid but approximate formula to yield the ionic thermal conductivity of warm dense Al with DPMD. Additionally, we recommend the use of a more complete formula for those systems whose ionic thermal conductivity is an important part of the total thermal conductivity.

Equation (7) can also be written asκI=∫0+∞CJ(t)dt,(9)from which the autocorrelation function of the heat current is defined asCJ(t)=13ΩkBT2〈Jq(t)⋅Jq(0)〉.(10)

We first perform FPMD simulations of Al based on KSDFT and OFDFT at temperatures of 0.5 eV, 1.0 eV, and 5.0 eV. The PBE XC functional is used, and the FPMD trajectory length is 10 ps. The cell contains 64 Al atoms. Two DP models, named DP-KS and DP-OF, are trained, based on the KSDFT and OFDFT trajectories, respectively. Note that the accuracy of the DP models in describing warm dense Al at the three temperatures has been demonstrated in our previous work,⁵⁹ where we showed that DPMD has excellent accuracy in reproducing structural and dynamical properties, including radial distribution functions and static and dynamic structure factors.

Here, we first focus on the frequency-dependent Onsager coefficients L_mn(ω). Figure 2 shows the L₁₁(ω), L₁₂(ω), and L₂₂(ω) computed via different methods at 0.5 eV; L₂₁(ω) is not shown, since L₂₁(ω) = L₁₂(ω), as derived from Eq. (4). We can see that for L₁₁(ω) and L₂₂(ω) all four methods yield very similar results, but for L₁₂(ω) there are deviations at low frequencies. The main reason for this is the greater sensitivity of L₁₂(ω) to the number of ionic configurations, and 20 snapshots are not enough to provide good convergence. We also run a test with 40 snapshots, and this improves the results, as can be seen in Fig. 2(d). However, the resulting electronic thermal conductivity depends mainly on L₂₂(ω) at the temperatures considered in this work, and so the small deviations in L₁₂(ω) do not affect the computed κ_e(ω), as will be shown next.

The results for the frequency-dependent electronic thermal conductivity κ_e(ω) computed using the four different methods are illustrated in Fig. 3. We can see that the OFDFT results agree well with the KSDFT ones, suggesting that OFDFT has the same accuracy as KSDFT in yielding atomic configurations for subsequent computations of κ_e(ω) via the KG method, even though some approximations are made with regard to the electron kinetic energy⁶⁵ and the local pseudopotential⁸⁵ within the framework of OFDFT. Impressively, the DP models trained from FPMD trajectories yield almost identical κ_e(ω) when compared with the DFT results, which proves that the DP models can yield highly accurate atomic configurations for subsequent calculations of κ_e(ω). In conclusion, the input atomic configurations for the calculation of κ_e(ω) can be generated by efficient DPMD models without loss of accuracy, which is beneficial for simulating a large number of atoms to mitigate size effects. Although the preparation of the training data in the DPMD model requires additional computational resources, we find that running FPMD with a cell consisting of 64 atoms is sufficient to generate reliable DPMD models. Therefore, additional computational costs are saved once larger numbers of atoms are adopted in the linear-scaling DPMD method.⁵⁹

Previous work^30,33,38 has shown that the electronic thermal conductivity κ_e depends strongly on both the number of k-points and the number of atoms. Additionally, the broadening parameter should be chosen appropriately to yield a meaningful κ_e. Here, we systematically investigate these issues by adopting the DPMD model to generate atomic configurations for subsequent calculations of κ_e. The DPMD model is trained from OFDFT-based MD trajectory with the PBE XC functional. By utilizing the snapshots from the DPMD trajectory, we obtain κ_e by using the KG formula. Furthermore, we study the effects of different XC functionals and pseudopotentials in affecting κ_e.

We first investigate the convergence of κ_e(ω) with respect to different numbers of k-points. Figure 4 illustrates an example of warm dense Al in a 64-atom cell at a temperature of 1.0 eV. The k-point meshes are chosen to be 1 × 1 × 1, 2 × 2 × 2, 3 × 3 × 3, and 4 × 4 × 4 in the calculations of κ_e(ω). The results show that κ_e(ω) converges when 3 × 3 × 3 k-points are used. In fact, the number of k-points required for convergence of κ_e(ω) varies with different numbers of atoms and at different temperatures. Table I lists the sizes of k-points that are needed to converge κ_e(ω) with numbers of atoms ranging from 16 to 256 at temperatures of 0.5 eV, 1.0 eV, and 5.0 eV. We find that the number of k-points required exhibits a trend to decrease at higher temperatures. For instance, the numbers of k-point samplings required for a 32-atom cell at temperatures of 0.5 eV, 1.0 eV, and 5.0 eV are 6 × 6 × 6, 5 × 5 × 5, and 3 × 3 × 3, respectively. This can be understood from the fact that the Fermi–Dirac distribution of electrons around the Fermi energy becomes sharper at low temperatures. Therefore, a dense mesh of k-points is needed to represent the detailed structures around the Fermi surface at low temperatures.

The value of κ_e(ω) converges not only for a sufficient number of k-points but also for a sufficient number of atoms in the simulation cell. To demonstrate this, we plot κ_e(ω) with respect to different numbers of atoms in Fig. 5. For each size of cell, the number of k-points is chosen large enough to converge κ_e(ω), as listed in Table I. We find the following. First, for temperatures of 0.5 eV, 1.0 eV, and 5.0 eV, we find that a 256-atom system is large enough to converge κ_e(ω). Second, in most cases, κ_e(ω) does not monotonically decrease but has a peak at low frequencies ranging from 0.3 eV to 0.6 eV. This peak moves toward ω = 0 when a larger number of atoms is used, and almost disappears in large cells such as the 432-atom cell at 0.5 eV. Similar to what was found in previous work,³⁸ this phenomenon is caused by size effects. Third, κ_e(ω) converges more rapidly with respect to the number of atoms at high temperatures, suggesting that size effects are less significant at higher temperatures. For example, the values of κ_e for the 16- and 64-atom systems at 0.5 eV are 59.7% (195.5 W m⁻¹ K⁻¹) and 15.6% (409.3 W m⁻¹ K⁻¹) lower than the value for the 256-atom system (485.1 W m⁻¹ K⁻¹). Meanwhile, the values of κ_e for the 16- and 64-atom systems at 5.0 eV are only 47.7% (669.7 W m⁻¹ K⁻¹) and 7.5% (1182.2 W m⁻¹ K⁻¹) lower than the value for a 256-atom system (1281.7 W m⁻¹ K⁻¹).

We perform a further analysis to elucidate the origin of the size effects in computations of κ_e(ω). As Eq. (4) shows, for a given energy interval ϵ_ik −ϵ_jk with electronic states i and j at a specific k point, κ_e(ω) should converge provided there are a sufficient number of electronic eigenstates. However, it can be seen from Fig. 5 that the computed κ_e(ω) becomes substantially larger at low frequencies as the number of atoms in the simulation cell increases. These results imply that even when there are enough k-points for convergence to be reached, if only a small number of atoms is used, then the small energy intervals obtained from the finite-size DFT calculations do not provide sufficient information to allow evaluation of κ_e(ω).

To clarify this issue, we define an energy interval distribution function (EIDF) asg(E)=1Np∑i>j,kW(k)δ(ϵik−ϵjk−E),(11)where W(k) represents the weights of the k-points used in the DFT calculations, ϵ_ik and ϵ_jk are the eigenvalues, and N_p is a normalization factor. We choose warm dense Al systems at temperatures of 0.5 eV and 5.0 eV, with the selected energy intervals computed from the bands occupied by 3s and 3p electrons, which can be identified from the density of states (DOS), as illustrated in Fig. 6. The energy intervals satisfy the condition that E < 1.0 eV. In addition, we consider the bands within 6.0 eV and 50.5 eV above the chemical potential μ for temperatures of 0.5 eV and 5 eV, respectively. The results for g(E) with respect to different numbers of atoms ranging from 16 to 432 are illustrated in Fig. 7. The EIDF g(E) becomes larger for small E as the number of atoms increases, and it converges when the number of atoms reaches 256 for both temperatures 0.5 eV and 5.0 eV. The increase in g(E) at small E implies that more energy eigenvalues with values close to one another appear in a larger cell with more atoms. In other words, because of the finite number of atoms used in the simulation cell, the energy levels are discretized to some extent, resulting in a lack of energy intervals, especially for those small values. In particular, these small energy intervals are of considerable importance in the evaluation of κ_e(ω) when ω → 0, as shown in Eq. (4).

The FWHM broadening parameter σ that appears in the δ(E) function in Eq. (6) substantially affects the resulting electronic thermal conductivity κ_e(ω) when ω → 0. We therefore investigate the influence of the choice of σ on the computed κ_e(ω) by analyzing the EIDF g(E). Taking Al at 0.5 eV as an example, we plot in Fig. 8 both g(E) and κ_e(ω) of a 256-atom cell, which is large enough to converge g(E) or κ_e(ω), as demonstrated above. When the broadening effect is small (σ = 0.01 eV), g(E) and κ_e(ω) decrease dramatically below 0.2 eV, which is a consequence of the discrete band energies. Thus, a suitable σ needs to be chosen to compensate for the discrete band energies. For instance, the value of g(E) at E = 0 increases from 0.22 eV⁻¹ (σ = 0.01 eV) to 0.99 eV⁻¹ (σ = 0.4 eV). However, g(E) becomes saturated if too large a σ is applied, resulting in overcorrection of κ_e(ω) such that the curve at frequencies lower than 0.6 eV decreases. Therefore, to compensate for the discrete band energies while avoiding overcorrection, we choose σ to be 0.4 eV for warm dense Al at all of the temperatures considered in this work. It is worth mentioning that a majority of studies, including this one, treat the broadening parameter as an adjustable variable to yield a better characterized curve of frequency-dependent thermal conductivity.^{18,29–31,33,34,38} Although some physical criteria have been proposed,⁸⁶ no definite conclusions have yet been drawn regarding the choice of the adjustable variable. In this work, we choose the adjustable parameter in terms of sufficiently small energy intervals, which are possible when a sufficiently large supercell is used in the DFT calculations. κ_e at zero frequency is obtained by linear extrapolation and is illustrated in Fig. 9.

We study the influences of the LDA and PBE XC functionals on the computed κ_e(ω) by first validating the atomic configurations generated by FPMD simulations. Specifically, atomic configurations are chosen from two 256-atom DPMD trajectories, which are generated by two DP models trained from OFDFT with the LDA and PBE XC functionals. We then adopt the KG method using the PP1 pseudopotentials generated with the same XC functional and obtain κ_e at temperatures of 0.5 eV, 1.0 eV, and 5.0 eV. As shown in Fig. 9 and listed in Table II, the κ_e values obtained from the PBE XC functional are 485.1 W m⁻¹ K⁻¹, 764.1 W m⁻¹ K⁻¹, and 1281.7 W m⁻¹ K⁻¹ at temperatures of 0.5 eV, 1.0 eV, and 5.0 eV, respectively, while the κ_e values from the LDA XC functional are 1.6% lower (477.3 W m⁻¹ K⁻¹), 1.3% higher (773.8 W m⁻¹ K⁻¹), and 2.8% lower (1246.4 W m⁻¹ K⁻¹) than those from the PBE XC functional at 0.5 eV, 1.0 eV, and 5.0 eV, respectively. Therefore, we conclude that the LDA XC functional yields almost the same κ_e values as the PBE XC functional. Note that in a recent work⁸⁷ that adopted the HSE XC functional, there were some differences in the electronic thermal and electrical conductivities of warm dense Al compared with those obtained from the PBE XC functional.

We investigate how norm-conserving pseudopotentials affect the computed electronic thermal conductivity κ_e. First of all, Fig. 6 shows the DOS of two types of pseudopotentials (PP1 and PP2) at temperatures of 0.5 eV and 5.0 eV, and we can see that the two pseudopotentials yield similar DOS of 3s3p electrons. Next, Fig. 9 illustrates the computed κ_e from two types of pseudopotentials and the computed data from Knyazev and Levashov,³¹ Vlcěk et al.,³⁴ and Witte et al.,⁸⁷ as well as the experimental data from McKelvey et al.⁸⁸ The values of κ_e from the two types of pseudopotentials are also listed in Table II.

We find the following. First, the DP-OF results agree reasonably well with the DP-KS ones, as has previously been shown in Fig. 3. For example, DP-KS with the PP1 pseudopotential and the PBE XC functional yields κ_e = 466.5 W m⁻¹ K⁻¹ at 0.5 eV, which is 3.8% lower than that from DP-OF (485.1 W m⁻¹ K⁻¹) at the same temperature, and the relative difference decreases to 1.9% at 5.0 eV. These results indicate that the OFDFT is appropriate for study of the electronic thermal conductivity of warm dense Al in the temperature range from 0.5 eV to 5.0 eV. Second, our calculations with the PP1 and PP2 pseudopotentials yield similar κ_e values of around 480 W m⁻¹ K⁻¹ and 770 W m⁻¹ K⁻¹ at temperatures of 0.5 eV and 1.0 eV, respectively. The results are consistent with those of Knyazev and Levashov, Vlcěk et al., and Witte et al. However, the κ_e values from the two pseudopotentials at a temperature of 5.0 eV deviate. For instance, the result from DP-OF (PBE) with the PP2 pseudopotential is 1604.6 W m⁻¹ K⁻¹, while the result with the PP1 pseudopotential is only 1281.7 W m⁻¹ K⁻¹, which is 20.1% lower than the former. The κ_e values from PP1 are close to those from Witte et al., while the values from PP2 are consistent with Knyazev and Levashov’s data. Note that Witte et al. used a PAW potential with 11 valence electrons and the PBE XC functional, whereas Knyazev and Levashov adopted an ultrasoft pseudopotential with three valence electrons and the LDA XC functional. It is also worth mentioning that Witte et al. used a 64-atom cell, whereas Knyazev and Levashov used a 256-atom cell, and so, according to our analysis, size effects may be present in the former work. In general, both values of electronic thermal conductivity lie within the range of the experimental data obtained by McKelvey et al. Even so, our results demonstrate that different pseudopotentials may substantially affect the results for κ_e at high temperatures. One possible reason for the deviation of κ_e at 5.0 eV is the number of electrons included in the pseudopotentials. However, it is also possible that the deviation comes from the fact that the nonlocal potential correction^77,78 in Eq. (4) is ignored. Therefore, nonlocal corrections must be considered when calculating the conductivity via the KG formula, and these will be the subject of our future work.

The ionic thermal conductivity of warm dense Al can be evaluated using the GK formula, since the atomic energies are available in the DPMD method. However, the computed ionic thermal conductivity may be affected by trajectory length and system size. We therefore study the convergence of the ionic thermal conductivity with respect to different lengths of trajectories and system sizes. We first test the convergence of the autocorrelation function C_J(t) in Eq. (10) with respect to different lengths of trajectories, and the results are shown in Fig. 10. A 10 648-atom Al system is considered, with four different lengths of trajectories, namely, 25 ps, 125 ps, 250 ps, and 500 ps, and the DPMD model is trained from OFDFT with the PBE XC functional at a temperature of 0.5 eV. As shown in Fig. 10, the C_J(t) curves abruptly decay within the first 0.1 ps and oscillate with respect to time t. We can see that the oscillations are strongly affected by the length of simulation time. For example, the C_J(t) curve obtained from the 25 ps trajectory exhibits substantially larger oscillations than the curves obtained from the other three trajectories, while the best convergence is achieved with the 250 ps trajectory. These results suggest that to obtain a converged C_J(t) for warm dense Al, a few hundreds of picoseconds are required, even for a system with more than 10⁴ atoms, which is beyond the capability of FPMD simulations but can be realized by DPMD simulations.

Next, to examine size effects on the ionic thermal conductivity, we run simulations for 500 ps for 12 different system sizes, with the number of atoms per cell ranging from 16 to 65 536. The results are shown in Fig. 11. Since C_J(t) cannot strictly reach zero, the correlation time t is truncated in the integral of C_J(t) in Eq. (9). In practice, we choose multiple truncations of t ranging from 0.5 ps to 1.5 ps and compute the error bars with the maximum and minimum integral values, which are also shown in Fig. 11. We find that the ionic thermal conductivity increases with increasing system size in the range from 16 to 1024 atoms per cell at all three temperatures considered for warm dense Al. With further increase in system size, the ionic thermal conductivity begins to oscillate until the largest system size is reached (65 536 atoms). We can conclude that at least a 1024-atom system should be adopted and that better convergence of ionic thermal conductivity can be achieved if a larger system is used. A 65 536-atom cell is adopted to compute the ionic thermal conductivity, and the results are shown in Table II. The ionic thermal conductivity of warm dense Al is around 1 W m⁻¹ K⁻¹–2 W m⁻¹ K⁻¹, which is more than two orders of magnitudes smaller than the electronic thermal conductivity. Additionally, it should be noted that the PBE and LDA XC functionals yield similar values for the ionic thermal conductivity.

We have proposed a method that combines DPMD and DFT to calculate both electronic and ionic thermal conductivities of materials, with the DP models being trained from DFT-based MD trajectories. The resulting DP models accurately reproduce the properties predicted by DFT. In addition, they can be utilized to efficiently simulate a large cell consisting of hundreds of atoms, which largely mitigates the size effects caused by periodic boundary conditions. By taking the atomic configurations from DPMD trajectories, we can use the eigenvalues and eigenstates of a given system obtained from DFT solutions and employ the KG formula to compute the electronic thermal conductivity. In addition, the DP models yield atomic energies, which are not available in the traditional DFT method. By using the atomic energies to evaluate the ionic thermal conductivity, both electronic and ionic contributions to the thermal conductivity can be obtained for a given material.

We took warm dense Al as an example and thoroughly studied its thermal conductivity. Expensive FPMD simulations of large systems can be replaced by DPMD simulations with much smaller computational resources. We first computed the temperature-dependent electronic thermal conductivities of warm dense Al from 0.5 eV to 5.0 eV at a density of 2.7 g/cm³ with snapshots from OFDFT, KSDFT, and DPMD, and the three methods yielded almost the same results, demonstrating that the DPMD method has similar accuracy to FPMD simulations. We then systematically investigated convergence issues with respect to the number of k-points, the number of atoms, the broadening parameter, the exchange-correlation functionals, and the pseudopotentials. A 256-atom system was found to be large enough to converge the electronic thermal conductivity. The broadening parameter was chosen to be 0.4 eV according to our analysis of the energy interval distribution function. We found that both LDA and PBE XC functionals yielded similar results for the electronic thermal conductivity. However, the choice of pseudopotentials may substantially affect the resulting electronic thermal conductivity. We also computed the ionic thermal conductivity with DPMD and the GK method and investigated convergence issues with respect to trajectory length and system size. We found that the ionic thermal conductivity of warm dense Al was much smaller than its electronic thermal conductivity. In summary, the DPMD method provides promising accuracy and efficiency in studying both electronic and ionic thermal conductivities of warm dense Al and should be considered for future work on modeling transport properties of WDM.