Thermal Transport by Electrons and Ions in Warm Dense Aluminum: A Combined Density Functional Theory and Deep Potential Study

We propose an efficient scheme, which combines density functional theory (DFT) with deep potentials (DP), to systematically study the convergence issues of the computed electronic thermal conductivity of warm dense Al (2.7 g/cm$^3$, temperatures ranging from 0.5 to 5.0 eV) with respect to the number of $k$-points, the number of atoms, the broadening parameter, the exchange-correlation functionals and the pseudopotentials. Furthermore, the ionic thermal conductivity is obtained by the Green-Kubo method in conjunction with DP molecular dynamics simulations, and we study the size effects in affecting the ionic thermal conductivity. This work demonstrates that the proposed method is efficient in evaluating both electronic and ionic thermal conductivities of materials.


I. INTRODUCTION
Warm dense matter (WDM) is a state of matter lying between condensed matter and plasma, which consists of strongly coupled ions and partially degenerated electrons.WDM exists in the interior of giant planets 1,2 or the crust of white dwarf and neutron stars 3,4 and can be generated through laboratory experiments such as diamond anvil cell, 5 gas gun 6,7 and high power lasers. 7,8WDM also plays an important role in the inertial confinement fusion (ICF). 9In this regard, it is crucial to understand properties of WDM such as the equation of state, optical and transport properties.7][18][19][20] Thermal conductivity is one of the upmost important properties of warm dense matter.2][23] It is a key parameter in the hydrodynamic instability growth on the National Ignition Facility capsule design. 24It also plays a key role in simulations of the interactions between laser and a metal target. 25Besides, the thermal conductivity largely affects predictions of ICF implosions in hydrosimulations. 26e thermal conductivity includes both electronic and ionic contributions.0][31][32][33][34][35] Typically, the κ e value is averaged over 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,37sides, the size effects may substantially affect κ e , as has been demonstrated in previous works. 30,33,38For computations of the ionic thermal conductivity κ I , we focus on utilizing the Green-Kubo (GK) formula, 27,39,40 where κ I is expressed by the heat flux auto-correlation function requiring the energy and the virial tensor of each atom from the molecular dynamics trajectory.5][46][47] Besides, a long molecular dynamics trajectory may be needed to evaluate κ I , which poses another challenge issue for DFT simulations.
While current experimental techniques only measure the total thermal conductivity, the first-principles methods have become an ideal tool to separately yield electronic and ionic contributions; however, only a few works that adopt first-principles methods have focused on transport properties by considering both contributions. 32,48,49For some metals, the ionic contribution is not important, 48 but this observation does not hold for all metals such as tungsten. 49Besides, the first-principles methods are computationally expensive especially when a large system or a long trajectory is considered.In this regard, the community awaits for an efficient and accurate method that can yield both electronic and ionic thermal conductivities of materials.Herein, with the above KG and GK formulas, we propose a combined DFT and deep potential (DP) method for the purpose and take warm dense Al as an example.The recently developed DP molecular dynamics (DPMD) [50][51][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 the traditional DFT, DPMD has a much higher efficiency while keeps the ab initio accuracy.In addition, the linear-scaling DPMD can be parallelized to simulate hundreds of millions of atoms. 53,54The DP method has been adopted in a variety of applications, such as crystallization of silicon, 55 high entropy materials, 56 isotope effects in liquid water, 57,58 and warm dense matter, 59,60 etc.
In this work, we demonstrate that the DP method in conjunction with DFT can be adopted to obtain both electronic and ionic thermal conductivities for 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, 1.0 and 5.0 eV with a density of 2.7 g/cm 3 .
The electronic thermal conductivity can be accurately computed via the KG method based on the DPMD trajectories as compared to the FPMD trajectories.Importantly, we systematically investigate the convergence issues such as the number of k-points, the number of atoms, the broadening parameter, the exchange-correlation functionals, and the pseudopotentials in affecting the electronic thermal conductivity with the aid of DPMD simulations.Furthermore, the ionic thermal conductivity can be obtained via DPMD and the convergence is studied with different sizes of systems and different lengths of trajectories.
The manuscript is organized as follows.In section II, we briefly introduce KSDFT, OFDFT, and DPMD and the setups of simulations.We also briefly introduce the Kubo-Greenwood formula and the Green-Kubo formula that used in this work.In section III, we first show the computed electronic thermal conductivity from the above three methods.
The convergence issues of the electronic thermal conductivity are then thoroughly discussed.
Finally, we show the results of the ionic thermal conductivity of warm dense Al.We conclude our works in Section IV.

A. Density Functional Theory
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 10 point out that the electron density which minimizes the total energy is the ground-state density while the energy is the ground-state energy. 10According to different treatments with the kinetic energy of electrons, two methods appear, i.e., KSDFT 11 and OFDFT. 12The kinetic energy of electron in KSDFT does not include the electron density explicitly but is evaluated from the ground-state wave functions of electrons obtained by the self-consistent iterations.On the other hand, OFDFT approximately denotes the kinetic energy of electrons as an explicit density functional, [61][62][63][64][65] which enables the direct search for the ground-state density and energy.As a result, KSDFT has a higher accuracy but OFDFT is several orders more efficient than KSDFT, especially for large systems.Besides, Mermin extends DFT to finite temperatures, 66 which has been widely used to describe electrons in WDM. 18,29 ran 64-atom Born-Oppenheimer molecular dynamics (BOMD) simulations with KSDFT by using the QUANTUM ESPRESSO 5.4 package. 67The Perdew-Burke-Ernzerhof (PBE) exchange-correlation (XC) functional 68 was used.Projector augmented-wave (PAW) potential 69,70 was adopted with three valence electrons and a cutoff radius of 1.38 Å.The plane wave cutoff energy was set to 20 Ry for temperatures of 0.5 and 1.0 eV and 30 Ry for 5.0 eV.We only used the gamma k-point.Periodic boundary conditions were used.Besides, the Andersen thermostat 71 was employed with the NVT ensemble and the trajectory length was 10 ps with a time step of 1.0 fs.
We also performed 108-atom BOMD simulations with OFDFT by utilizing the PROFESS 3.0 package. 72Both PBE 68 and LDA 11 XC functionals were used, in together with the Wang-Teter (WT) KEDF 65 and the Nosé-Hoover thermostat 73,74 in the NVT ensemble.We ran OFDFT simulations for 10 ps with a time step of 0.25 fs and the energy cutoff was set to 900 eV.Periodic boundary conditions were utilized.

B. Deep Potential Molecular Dynamics
The DP method [50][51][52] learns the dependence of the total energy on the coordinates of atoms in a system and build a deep neural network (DNN) model, which predicts the potential energy and force of each atom.In the training process, the total potential energy E tot is decomposed into energies of different atoms: where E i is the potential energy of atom i.In general, for each atom i, the mapping is established through DNN between E i and the atomic coordinates of its neighboring atoms within a cutoff radius r c .Specifically, the DNN model consists of the embedding network and the fitting network. 50The embedding network imposes constrains to atoms, enabling atomic coordinates to obey the translation, rotation and permutation symmetries.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 the loss function.The resulting DNN-based model can be used to simulate a large system with the efficiency comparable to empirical force fields.
In this work, we adopted the DNN-based models trained from either KSDFT or OFDFT trajectories with the DeePMD-kit package. 51With the purpose to study the size effects of electronic thermal conductivity, a series of cubic cells consisting of 16, 32, 64, 108, 216 and 256 atoms were simulated for 10 ps with a time step of 0.25 fs.Besides, we ran systems of 16, 32, 64, 108, 256, 1024, 5488, 8192, 10648, 16384, 32000 and 65536 atoms in order to investigate the convergence of ionic thermal conductivity.The length of these trajectories is 500 ps with a time step of 0.25 fs.We employed the Nosé-Hoover thermostat 73,74 in the NVT ensemble.Periodic boundary conditions were adopted.All of the DPMD simulations were performed by the modified LAMMPS package. 75

C. Kubo-Greenwood Formula
Electronic thermal conductivity κ e is calculated by the Onsager coefficients L mn as where T is temperature and e is the charge of electrons, and L mn is obtained by the frequencydependent Onsager coefficients L mn (ω) as We follow the Kubo-Greenwood formula derived by Holst et al., 76 and L mn (ω) takes the form of where m e is the mass of electrons, Ω is the volume of cell, W (k) represents the weight of k points in the Brillouin zone, µ is the chemical potential, Ψ i,k represents the wave function of the i-th band with the eigenvalue i,k and f being the Fermi-Dirac function.
We used the chemical potential µ instead of the enthalpy per atom, which does not affect the results of electronic thermal conductivity in one-component systems. 76[79] However, the non-local errors were demonstrated to be sufficiently small (about 3.3%) for liquid Al, 78 which are considered to be reasonably small in this work.In future works, we will include the additional term caused by non-local pseudopotentials.We also define the frequency-dependent electronic thermal conductivity as In practical usage of the KG method, the delta function in Eq. ( 4) needs to be broadened.
We adopte a Gaussian function 29 and the delta function takes the form of Here ∆E controls the full width at half maximum (FWHM, denoted as σ) of Gaussian function with the relation of σ ≈ 2.3548∆E.
The KG method needs computed eigenvalues and wave functions from DFT solutions of given atomic configurations.In practice, we selected 5-20 atomic configurations from the last 2-ps MD trajectories with a time interval of 0.1 ps.We used both PBE and LDA XC functionals and the associated norm conserving (NC) pseudopotentials in order to test the influences of XC functionals and pseudopotentials on the resulting electronic thermal conductivity.We adopted two NC pseudopotentials for Al, which are referred as PP1 and PP2.The PP1 pseudopotential was generated with the optimized normconserving Vanderbilt pseudopotential method via the ONCVPSP package. 80,81We used 11 valence electrons and a cutoff radius of 0.50 Å. Calculations of κ e in a 64-atom cell involved 1770, 2100 and 5625 bands at temperatures of 0.5, 1.0 and 5.0 eV, respectively.The PP2 pseudopotential was generated through the PSlibrary package. 82We used the Troullier-Martins method 83 and the cutoff radius was set to 1.38 Å.We chose 3 valence electrons for each atom.We selected 720, 1100 and 4800 bands for calculations of κ e at temperatures of 0.5, 1.0 and 5.0 eV, respectively.The plane wave cutoff energies of both pseudopotentials were set to 50 Ry.Generally, the PP1 pseudopotential was used in most cases and the PP2 pseudopotential was only used to compare the effects of different pseudopotentials on κ e .Fig. 1 shows the computed averaged κ e (ω) with respect to different numbers of atomic configurations from 108-atom DPMD trajectories at 0.5 and 5.0 eV.The DP model was 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 (ω) is easier to converge at the relatively lower temperature of 0.5 eV.Therefore, we chose 20 atomic configurations for cells with 108 atoms or less at all of the three temperatures.For systems with the number of atoms larger than 108, we respectively selected 5 and 10 atomic configurations for temperatures smaller than 5.0 eV and equal to 5.0 eV unless otherwise specified.

D. Green-Kubo Formula
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 40 with the formula of where Ω is the volume of cell, T is temperature, k B is the Boltzmann constant, • • • is 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 which takes the form of Here v i is the velocity of the i-th atom, while ε i is the energy of the i-th atom including ionic kinetic energy and potential energy, where the potential energy is directly obtained by DNN.F ij is the force acting on the i-th atom due to the presence of the j-th atom with r ij defined as r i − r j .It is worth mentioning that Eq. ( 8) is valid for 2-body potentials but yields about 20% error for J q when a many-body potential (e.g. the deep potential) is adopted. 84Importantly, 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 counterpart.Therefore, we still consider Eq. ( 8) as a valid but approximated formula to yield the ionic thermal conductivity of warm dense Al with DPMD.Additionally, we recommend the usage of a more complete formula for those systems whose ionic thermal conductivity is an important part of the total thermal conductivity.
Eq. ( 7) can also be written as from which the auto-correlation function of heat current is defined as III. RESULTS

A. Accuracy of DP models
We first performed FPMD simulations of Al based on KSDFT and OFDFT at temperatures of 0.5, 1.0 and 5.0 eV.The PBE XC functional was used and the FPMD trajectory length was 10 ps.The cell contained 64 Al atoms.Two DP models named DP-KS and DP-OF were 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 have been demonstrated in our previous work, 59 where we show that DPMD has an excellent accuracy in reproducing structural and dynamical properties including radial distribution functions, static structure factors, and dynamic structure factors.
Here, we first focus on the frequency-dependent Onsager coefficients L mn (ω).Fig. 2 shows the computed L 11 (ω), L 12 (ω), and L 22 (ω) via different methods at 0.5 eV; L 21 (ω) is not illustrated since L 21 (ω) = L 12 (ω) as derived from Eq. ( 4).We can see that all of the four methods yield very close L 11 (ω) and L 22 (ω) except for L 12 (ω), where slightly larger differences are found at low frequencies.The main reason is that L 12 (ω) is more sensitive to the number of ionic configurations and 20 snapshots are not enough to converge it well.
We also test 40 snapshots and the results are improved, as illustrated in Fig. 2(d).However, the resulting electronic thermal conductivity mainly depends on L 22 (ω) at temperatures considered in this work and the small differences of L 12 (ω) do not affect the computed κ e (ω), as will be shown next.
The computed frequency-dependent electronic thermal conductivity κ e (ω) utilizing the abovementioned four different methods are illustrated in Fig. 3.We can see that the OFDFT results agree well with the KSDFT ones, suggesting OFDFT has the same accuracy as KSDFT to yield atomic configurations for subsequent computations of κ e (ω) using the KG method, even though there are some approximations on the kinetic energy of electrons 65 and the local pseudopotential 85 within the framework of OFDFT.Impressively, the DP models trained from FPMD trajectories yield almost identical κ e (ω) when compared to 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 losing 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 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. 59

B. Convergence of Electronic Thermal Conductivity
Previous works 30,33,38 have 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 properly chosen in order to yield a meaningful κ e .Herein, we systematically investigate the above issues by adopting the DPMD model to generate atomic configurations for subsequent calculations of κ e .The DPMD model was 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 .

Number of k-points
We first investigate the convergence of κ e (ω) with respect to different k-points.Fig. 4 illustrates an example of warm dense Al in a 64-atom cell at a temperature of 1.0 eV.The k-point meshes were chosen to be 1

Number of atoms
The value of κ e (ω) converges not only with enough number of k-points but also with sufficient number of atoms in the simulation cell.To demonstrate this point, 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 have the following findings.First, for the temperatures of 0.5, 1.0 and 5.0 eV, we find that a 256atom system is large enough to converge κ e (ω).Second, most of κ e (ω) do not monotonically decrease but have peaks at low frequencies ranging from 0.3 to 0.6 eV.These peaks move towards ω=0 when a larger number of atoms are utilized, and almost disappear in large cells such as the 432-atom cell at 0.5 eV.Similar to previous works, 38 this phenomenon is caused by the size effects.Third, κ e (ω) converges faster with respect to the number of atoms at high temperatures, suggesting that the size effects are less significant at higher temperatures.For example, κ e of the 16-and 64-atom systems at 0.5 eV are 59.7% (195.5 Wm −1 K −1 ) and 15.6% (409.3Wm −1 K −1 ) lower than the value from a 256-atom system (485.1 Wm −1 K −1 ).Meanwhile, κ e of the 16-and 64-atom systems at 5.0 eV are only 47.7% (669.7 Wm −1 K −1 ) and 7.5% (1182.2Wm −1 K −1 ) lower than the value from a 256-atom system (1281.7 Wm −1 K −1 ).
We perform further analysis to elucidate the origin of size effects in computations of κ e (ω), which is due to insufficient small energy intervals caused by limited sizes of systems.
As Eq. ( 4) shows, for a given energy interval ik − jk with electronic states of i and j at a specific k point, κ e (ω) should converge with enough electronic eigenstates; however, Fig. 5 illustrates that the computed κ e (ω) substantially becomes larger at low frequencies with increased number of atoms in the simulation cell.The results imply that the information of small energy intervals gained from finite-size DFT calculations are insufficient to evaluate κ e (ω) even when the number of k-points reaches convergence but a small number of atoms is adopted.To clarify this issue, we define an energy interval distribution function (EIDF) as where W (k) represents the weights of k-points used in DFT calculations, ik and jk are the eigenvalues and N p is a normalization factor.We chose warm dense Al systems at temperatures of 0.5 and 5.0 eV with the selected energy intervals computed from the bands occupied by 3s and 3p electrons, which can be identified from 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 and 50.5 eV above the chemical potential µ for temperatures of 0.5 and 5 eV, respectively.The results of g(E) with respect to different numbers of atoms ranging from 16 to 432 atoms are illustrated in Fig. 7.The EIDF g(E) becomes larger for small E as the number of atoms increases and converges when the number of atoms reaches 256 for both cases at 0.5 and 5.0 eV.The increase of g(E) at small E implies that more energy eigenvalues that have close values 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 the lack of energy intervals especially for those small values.Particularly, these small energy intervals are of significant importance in evaluating κ e (ω) when ω → 0 as shown in Eq. ( 4).

Broadening Parameter
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 choices of σ in influencing 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 caused by the discrete band energies.Thus, a suitable σ needs to be chosen to compensate the discrete band energies.For instance, the value of g(E) at E=0 increases from 0.22 eV −1 (σ=0.01 eV) to 0.99 eV −1 (σ=0.4 eV).However, g(E) becomes saturated if a too large σ is applied, resulting in overcorrected κ e (ω) that the curve at frequencies lower than 0.6 eV decreases.Therefore, in order to compensate the discrete band energies and avoid overcorrection at the same time, we choose σ to be 0.4 eV for warm dense Al at all of the temperatures considered in this work.34,38 Even though some physical criteria have been proposed, 86 the choice of the adjustable variable is still inconclusive.In this work, we choose the adjustable parameter in terms of sufficient small energy intervals, which are possible when a sufficiently large supercell is used in DFT calculations.The κ e at zero frequency is obtained by the linear extrapolation and illustrated in Fig. 9.

Exchange-correlation functionals
We study the influences of LDA and PBE XC functionals on the computed κ e (ω) by first validating the atomic configurations generated by FPMD simulations.Specifically, atomic configurations were chosen from two 256-atom DPMD trajectories, which were generated by two DP models trained from OFDFT with the LDA and PBE XC functionals.We then adopted the KG method by using the PP1 pseudopotentials generated with the same XC functional, and yielded κ e at temperatures of 0.5, 1.0, 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,764.1 and 1281.7 Wm −1 K −1 at temperatures of 0.5, 1.0 and 5.0 eV, respectively, while the κ e values from the LDA XC functional are 1.6% lower (477.3Wm −1 K −1 ), 1.3% higher (773.8 Wm −1 K −1 ) and 2.8% lower (1246.4Wm −1 K −1 ) than those of PBE at 0.5, 1.0 and 5.0 eV respectively.Therefore, we conclude that the LDA XC functional yields almost the same κ e values as PBE.Note that a recent work 87 adopted the HSE XC functional and showed some differences in electronic thermal and electrical conductivities of warm dense Al as compared with those from the PBE XC functional.

Pseduopotentials
We investigate how norm-conserving pseudopotentials affect the computed electronic thermal conductivity κ e .First of all, Fig. 6 shows DOS of two types of pseudopotentials (PP1 and PP2) at temperatures of 0.5 and 5.0 eV, where we see that the two pseudopotentials yield similar DOS of 3s3p electrons.Next, Fig. 9 illustrates the computed κ e from two types of pseudopotentials and those computed data from Knyazev et al., 31 Vlcěk et al. 34  We have the following findings.First, the DP-OF results agree reasonably well with the DP-KS ones, as have been previously shown in Fig. 3.For example, DP-KS with the PP1 pseudopotential and the PBE XC functional yields κ e =466.5 Wm −1 K −1 at 0.5 eV, which is 3.8% lower than that of DP-OF (485.1 Wm −1 K −1 ) at the same temperature and the relative difference decreases to 1.9% at 5.0 eV.The above results imply that the OFDFT is suitable to study the electronic thermal conductivity of warm dense Al ranging from 0.5 is the number of electrons included in the pseudopotentials.However, it is also possible that the deviation comes from the fact that the non-local potential correction 77,78 in Eq. ( 4) is ignored.Therefore, non-local corrections have to be considered when calculating the conductivity via the Kubo-Greenwood formula which is subject of future work.

C. Ionic Thermal Conductivity
The ionic thermal conductivity of warm dense Al can be evaluated by 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.In this regard, we 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 auto-correlation function C J (t) in Eq. ( 10) with respect to different lengths of trajectories and the results are shown in Fig. 10.A 10648-atom Al system was adopted with four different lengths of trajectories, i.e., 25, 125, 250, and 500 ps, and the DPMD model was trained from OFDFT with the PBE XC functional at a temperature of 0.5 eV.As illustrated in Fig. 10, the C J (t) curves abruptly decay within the first 0.1 ps and oscillate with respect to time t.We notice that the oscillations are largely affected by the length of simulation time.For example, C J (t) obtained from the 25-ps trajectory exhibits substantially larger oscillations than the other three trajectories, and the convergence is better achieved when the 250-ps trajectory is adopted.The above results suggest that in order to yield converged C J (t) for warm dense Al, a few hundreds of ps are required even for a system with more than ten thousand atoms, which is beyond the capability of FPMD simulations but can be realized by DPMD simulations.
Next, we run 12 different sizes of cells ranging from 16 to 65536 atoms for 500 ps to check the size effects on ionic thermal conductivity, and the results are illustrated in Fig. 11.
Since C J (t) cannot strictly reach zero, a truncation of the correlation time t is applied to the integration of C J (t) in Eq. ( 9).In practice, we choose multiple truncations of t ranging from 0.5 to 1.5 ps and computed 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 larger system sizes ranging from 16-to 1024 atoms at all of the three temperatures considered for warm dense Al.Next, the ionic thermal conductivity begins to oscillate until the largest system size adopted (65536 atoms).In this regard, we conclude that at least a 1024-atom system should be adopted and a better converged ionic thermal conductivity can be obtained if a larger size of system is used.A 65536-atom cell is utilized to compute the ionic thermal conductivity and the results are in Table II.The ionic thermal conductivity of warm dense Al is around 1-2 Wm −1 K −1 , which is more than two orders of magnitudes smaller than the electronic thermal conductivity counterpart.Additionally, both PBE and LDA XC functionals yield similar values for the ionic thermal conductivity.

IV. CONCLUSIONS
We propose a method that combines DPMD and DFT to calculate both electronic and ionic thermal conductivities of materials, and the DP models are trained from DFT-based MD trajectories.The resulting DP models accurately reproduce the properties as compared to those from DFT.In addition, the DP models can be utilized to efficiently simulate a large cell consisting of hundreds of atoms, which largely mitigate the size effects caused by periodic boundary conditions.Next, by using the atomic configurations from DPMD trajectories, one can used the eigenvalues and eigenstates of a given system obtained from DFT solutions, and adopt the Kubo-Greenwood 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 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 to 5.0 eV at a density of 2.7 g/cm 3 with snapshots from OFDFT, KSDFT and DPMD, and the three methods yielded almost the same results, demonstrating that the DPMD method owns similar accuracy as FPMD simulations.We then systematically investigated the convergence issues with respect to the number of k-points, the number of atoms, the broadening parameter, the exchangecorrelation 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 both LDA and PBE XC functionals yielded similar results for the electronic thermal conductivity.However, the choices of pseudopotentials may substantially affect the resulting electronic thermal conductivity.Furthermore, we also computed the ionic thermal conductivity with DPMD and the GK method, and investigated the convergence issues with respect to trajectory length and system size.We found the ionic thermal conductivity of warm dense Al is much smaller than its electronic thermal conductivity.In summary, the DPMD method provides a promising accuracy and efficiency in studying both electronic and ionic thermal conductivity of warm dense Al and should be considered for future work on modeling transport properties of WDM.

and 4 × 4 ×
4 in calculationsof κ e (ω).The results show that κ e (ω) converges when 3 × 3 × 3 k-points are used.In fact, the required number of k-points for convergence of κ e (ω) varies with different number of atoms and different temperatures.TableIlists 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, 1.0 and 5.0 eV.We find that the needed number of k-points exhibits a trend to decrease at higher temperatures.For instance, the needed k-point samplings of a 32-atom cell at temperatures of 0.5, 1.0 and 5.0 eV are 6 × 6 × 6, 5 × 5 × 5 and 3 × 3 × 3, respectively.This can be understood by the fact that the Fermi-Dirac distribution of electrons around the Fermi energy becomes more sharp at low temperatures.Therefore, a dense mesh of k-points is needed to represent the detailed structures around the Fermi surface at low temperatures.

to 5 .
0 eV.Second, our calculations with the PP1 and PP2 pseudopotentials yield similar κ e values of around 480 and 770 Wm −1 K −1 at 0.5 and 1.0 eV, respectively.The results are consistent with those from Knyazev et al., Vlcěk et al. and Witte et al.However, the κ e values from the two pseudopotentials at 5.0 eV deviate.For instance, the result of DP-OF (PBE) with the PP2 pseudopotential is 1604.6Wm −1 K −1 while the result with the PP1 pseudopotential is only 1281.7 Wm −1 K −1 , which is 20.1% lower than the former one.The κ e values from PP1 are close to those from Witte et al., while the κ e values from PP2 are consistent with the Knyazev et al. data.Note that Witte et al. utilized a PAW potential with 11 valence electrons and the PBE XC functional, while Knyazev et al. adopted an ultrasoft pseudopotential with 3 valence electrons and the LDA XC functional.It is also worth mentioning that a 64-atom cell is utilized in the work by Witte et al. while a 256-atom cell is used in the work by Knyazev et al.In this regard, the size effects may exist in the previous one according to our analysis.In general, both values of electronic thermal conductivity lie within the experimental data of McKelvey et al.Even so, our results demonstrate that different pseudopotentials may substantially affect the results of κ e at high temperatures.One possible reason that causes the deviation of κ e at 5.0 eV
and Witte et al., 87 as well as the experimental data from McKelvey et al.. 88 Besides, κ e from two types of pseudopotentials are also listed in Table.II.

TABLE II .
Electronic thermal conductivity κ e (Wm −1 K −1 ) and ionic thermal conductivity κ I (Wm −1 K −1 ) at temperatures T of 0.5, 1.0 and 5.0 eV.The results are computed from the DP-KS and DP-OF molecular dynamics trajectories.DP-KS and DP-OF refer to the DP models trained from KSDFT and OFDFT molecular dynamics trajectories, respectively.