Combining stochastic density functional theory with deep potential molecular dynamics to study warm dense matter

In traditional finite-temperature Kohn-Sham density functional theory (KSDFT), the well-known orbitals wall restricts the use of first-principles molecular dynamics methods at extremely high temperatures. However, stochastic density functional theory (SDFT) can overcome the limitation. Recently, SDFT and its related mixed stochastic-deterministic density functional theory, based on the plane-wave basis set, have been implemented in the first-principles electronic structure software ABACUS [Phys. Rev. B 106, 125132(2022)]. In this study, we combine SDFT with the Born-Oppenheimer molecular dynamics (BOMD) method to investigate systems with temperatures ranging from a few tens of eV to 1000 eV. Importantly, we train machine-learning-based interatomic models using the SDFT data and employ these deep potential models to simulate large-scale systems with long trajectories. Consequently, we compute and analyze the structural properties, dynamic properties, and transport coefficients of warm dense matter. The abovementioned methods offer a new approach with first-principles accuracy to tackle various properties of warm dense matter.

In traditional finite-temperature Kohn-Sham density functional theory (KSDFT), the partial occupation of a large number of high-energy KS eigenstates restricts the use of first-principles molecular dynamics methods at extremely high temperatures.However, stochastic density functional theory (SDFT) can overcome the limitation.Recently, SDFT and its related mixed stochastic-deterministic density functional theory, based on the plane-wave basis set, have been implemented in the first-principles electronic structure software ABACUS [Phys.Rev. B 106, 125132 (2022)].In this study, we combine SDFT with the Born-Oppenheimer molecular dynamics (BOMD) method to investigate systems with temperatures ranging from a few tens of eV to 1000 eV.Importantly, we train machine-learning-based interatomic models using the SDFT data and employ these deep potential models to simulate large-scale systems with long trajectories.Consequently, we compute and analyze the structural properties, dynamic properties, and transport coefficients of warm dense matter.

I. INTRODUCTION
Understanding materials in extremely hightemperature conditions, such as warm dense matter (WDM) and hot dense plasma (HDP), is essential for comprehending various objects, such as planetary objects 1 , inertial confinement fusion 2 , and high-intensity, high-energy laser pulse experiments 3 .However, the extremely high temperatures of WDM and HDP pose significant challenges in terms of measuring their properties experimentally or predicting them theoretically 4 .After decades of combined efforts from both experimental 5 and theoretical 4 perspectives, it has been recognized that an adequate quantum-mechanical description of electrons is essential for theoretical calculations to have the prediction power.In fact, various first-principles approaches based on different levels of approximations are available to address this issue.To name a few, the Kohn-Sham density functional theory (KSDFT) 6,7 , the path-integral Monte Carlo (PIMC) [8][9][10][11] , the orbital-free density functional theory (OFDFT) [12][13][14][15] , the extended first-principles molecular dynamics (Ext-FPMD) [16][17][18][19] , and the stochastic density functional theory (SDFT) [20][21][22][23][24] methods are most popularly employed.
KSDFT with the Mermin finite-temperature densityfunctional approach 25 is typically employed to compute properties of materials at low temperatures.However, when temperatures are elevated to the WDM regime, the partial occupation of a large number of high-energy KS eigenstates becomes a severe hurdle as the number of electronic states that need to be considered is proportional to T been successfully applied to study the equation of states for low Z elements 30 , PIMC is largely limited at lower temperatures 31 .Compared with KSDFT, the OFDFT method avoids diagonalizing the wave functions and owns a higher efficiency [12][13][14][15] .However, the applications of OFDFT are still limited by the lack of satisfactory accuracy in the kinetic energy density functional 32 .
The Ext-FPMD method has been proposed to evaluate systems at extremely high temperatures with the first-principles accuracy [16][17][18][19] .This method treats the wave functions of high-energy electrons as plane waves analytically, thereby avoiding the partial occupation of a large number of high-energy KS eigenstates.However, it is still challenging to adopt the Ext-FPMD method to investigate the electrical conductivity and electronic thermal conductivity of materials due to the lack of orbital information for high-energy electrons.Similarly, the SDFT scheme adopts stochastic orbitals in combination with the Chebyshev trace method to overcome the partial occupation of a large number of high-energy KS eigenstates limit.Moreover, the statistical errors can be effectively reduced when more stochastic orbitals are included or larger systems are utilized.In addition, the mixed stochastic-deterministic DFT (MDFT) combines the KS and stochastic orbitals and speeds up calculations 33 .Note that the SDFT and MDFT methods, based on the plane-wave basis set, can be used with the k-point sampling method and have recently been implemented in the first-principles package ABACUS [34][35][36] .Since both SDFT and MDFT employ stochastic wave functions, we will refer to these two methods as SDFT throughout the remainder of this paper.
Recently, machine-learning-assisted atomistic simulation methods have achieved remarkable success and attracted considerable attention [37][38][39][40][41][42][43][44][45][46][47][48] .In particular, deep neural networks are often adopted to learn the potential energy surfaces, which are formed by the relations between atomic configurations and their corresponding energies, forces, and stresses.Essentially, these neural- network-based models demonstrate high efficiency while maintaining ab initio accuracy, as they effectively bypass the need to solve quantum mechanics-based equations.Among them, the deep potential molecular dynamics (DPMD) method [49][50][51][52] achieves high performances 53,54 and stands out as a promising method to study WDM.For example, the DPMD method has been applied to study the structural and dynamic properties of warm dense aluminum 15,55 , the ion-ion dynamic structure factor of warm dense aluminum 55 , the thermal transport by electrons and ions of warm dense aluminum 56 , and the principal Hugoniot curves of warm dense beryllium 51 , etc.
To summarize, two main challenges exist for largescale first-principles simulations of WDM.First, the existence of the partial occupation of a large number of highenergy KS eigenstates results in computationally expensive simulations of WDM at high temperatures.Second, obtaining converged results for certain physical properties of WDM with a small number of atoms is difficult unless a large cell with a long MD trajectory is obtained.In this work, we first validate the accuracy of SDFT by analyzing the statistical errors from the stochastic orbitals and compare with results from the traditional KS-DFT method.Here, we select warm dense boron (B) and carbon (C) as benchmark systems.Second, we couple the Born-Oppenheimer molecular dynamics (BOMD) method with SDFT to simulate warm dense B. Specifically, we generate two DP models to describe B at a density of 2.46 g/cm 3 with two different temperatures (86 and 350 eV); the training data are obtained from SDFTbased BOMD simulations.Third, by performing DPMD simulations, we significantly extend the time and spatial scales of warm dense B and obtain converged data for certain physical properties.The above workflow is shown in Fig. 1.Our work demonstrates that combining SDFT with the deep potential (DP) method offers a promising route to simulate WDM over a wide range of temperatures.
The rest of the paper is organized as follows.Section II describes the computational methods utilized in this work.In Section III, we discuss the results obtained from SDFT and the DP model.Finally, we summarize our results in Section IV.

A. Stochastic and Mixed Stochastic-Deterministic Density Functional Theory
In the KSDFT framework 6,7 , the single-particle DFT Hamiltonian is defined as where the first term depicts the kinetic operator, VHartree denotes the Hartree potential, Vxc is the exchangecorrelation (XC) potential, and Vext is the potential of interactions between the electrons and the nuclei as well as other external fields.Within the finite-temperature Mermin-Kohn-Sham theory 25 , the electron density is given by where r is the position operator of the electron, the prefactor 2 accounts for the electron spin.The parameter ω k is the weight of the k point with N k being the number of k points, while N KS is the number of occupied orbitals with i being the index.The Fermi-Dirac distribution function takes the form of where µ is the chemical potential.Here, ϕ ik (r) and ε ik respectively depict the eigenfunction and eigenvalue of the self-consistent KS Hamiltonian In practice, solving ϕ ik (r) and ε ik is costly at high temperatures since the needed number of KS wave functions is proportional to T 3 2 21 .Given any orthogonal and complete basis set {ψ j }, a stochastic orbital χ nk in SDFT 20,21 can be defined as which satisfies where θ nk j is randomly generated by a uniform distribution between 0 and 1, and N sto is the number of stochastic orbitals.
In MDFT 33 , both deterministic orbitals ϕ ik and stochastic orbitals χnk are used, Here, the stochastic orbitals are defined to be orthogonal to the deterministic orbitals.The number of deterministic orbitals is set to N KS , which is typically chosen to be a subset of occupied states.In addition, both sets of orbitals satisfy the relation Then, the electron density is given by where f ĤDFT , µ, T is calculated by the Chebychev expansion 57 .If N sto = 0 or N KS = 0, the MDFT method changes to the standard KSDFT or SDFT methods, respectively.
Notably, using stochastic orbitals in practical calculations results in statistical errors since these orbitals only form a complete basis when the number of stochastic orbitals approaches infinity.As reported in previous studies 23,36,58 , the error caused by the stochastic orbitals is proportional to 1/ √ N sto .However, when periodic boundary conditions (PBCs) with the k-point sampling are considered and each k-point has N sto stochastic orbitals, the resulting σ s is proportional to 1/ √ N k N sto , suggesting that more k-points can reduce the stochastic errors.In order to evaluate the accuracy of SDFT, we perform both SDFT and KSDFT calculations for B and C bulk systems in extreme conditions.
In the first-principles MD simulations of WDM, the KSDFT couples with the dynamics of ions, usually through the BOMD method.Since the motions of ions are treated classically in the BOMD method, we need to evaluate the force of an atom i in the form of Here E is the sum of electrons' energy and ion-ion repulsion energy, and R i is the position of atom i.By utilizing the plane-wave basis set and norm-conserving pseudopotentials within the traditional KSDFT and SDFT methods, the force of an atom can be decomposed into three parts where F L i is the local pseudopotential force term, F N L i depicts the nonlocal pseudopotential force term, and F II i is the Ewald force term origins from the ion-ion interactions.Furthermore, stress is defined as where ϵ αβ is the strain with the spatial coordinates α and β.The kinetic energy term of electrons is σ T αβ .σ L αβ and σ N L αβ are the local and nonlocal pseudopotential terms, respectively.σ Hartree αβ is the Hartree and σ xc αβ is the exchange-correlation term.The Ewald term is σ II αβ .One can refer to Ref. 36 by Liu et al. for more details on implementing the total energy, the total free energy, forces, and stresses within the framework of SDFT in ABACUS.Note that the SDFT method employed in this work still uses plane wave basis, so it is still computationally demanding for a system consisting of a few hundred atoms.

B. Deep Potential Molecular Dynamics
In the DPMD method 50,59 , the total energy E of a system is expressed as a sum of atomic contributions, i.e., E = i E i , where the energy E i from atom i depends on an environment matrix R i , which includes the information of neighboring atoms of atom i within a cutoff radius.The DP model maps R i via an embedding neural network to a symmetry-preserving descriptor, and then the descriptor is mapped to a fitting neural network to yield E i .A loss function is utilized to optimize the parameters in the embedding and fitting networks to generate DP models.The loss function is defined as where N is the number of atoms, ϵ = E/N is the energy per atom, F i is the force acting on atom i, ξ is the virial tensor per atom, and ∆ denotes the difference between the training data and the results predicted by the DP model.In addition, p ϵ , p f , and p ξ are tunable prefactors.The stochastic gradient descent scheme Adam 60 is adopted to train the DP model.We separately train the data from SDFT-based BOMD trajectories of B at 86 and 350 eV, where the Perdew-Burke-Ernzerhof (PBE) 61 functional is used.As a result, we obtain two DP models of warm dense B at the two temperatures.To better characterize the warm dense B at 350 eV, we adopt the temperature-dependent deep potential (TDDP) method 51 to train the DP model.Note that the TDDP method, as shown in Fig. 1(b), introduces the electron temperature of the system into the fitting net, which is more suitable for high-temperature systems.
Both embedding and fitting neural networks contain three layers with the specific number of neurons being (25, 50, 100) and (120, 120, 120), respectively.The cutoff radius for each atom is chosen to be 6.0 Å.The inverse distance 1/r decays smoothly from 0.5 to 6.0 Å in order to remove the discontinuity introduced by the cutoff.Both DP models undergo training for 500,000 steps.Throughout the training process, the values of p ϵ , p f , and p ξ are gradually adjusted from 0.02 to 1, 1000 to 1, and 0.02 to 1, respectively.We also employ the DP compress technique to both DP models to accelerate the DPMD simulations, as described in the literature 62 .

A. Statistical Errors of SDFT
To analyze the statistical errors that arise from the SDFT method itself, we choose a 32-atom B system with a density of 2.46 g/cm 3 at the temperature of 350 eV.In addition, we employ the PBE 61 exchange-correlation functional.We note that at the temperature of 86 eV and higher temperatures, the pseudopotential of B is generated by the ONCVPSP 63 method with all of its 5 electrons.The cutoff radius for the pseudopotential is set to 0.7 Bohr to avoid overlaps of electron orbitals at high temperatures.In addition, we select an energy cutoff of For each data point, the RMSE is obtained via 9 independent SDFT calculations with different sets of stochastic orbitals.
180, 240, and 300 Ry for temperatures of 86, 350, and 1000 eV, respectively.Figure 2 shows the RMSE of SDFT verse the number of stochastic orbitals (N sto ) and the number of k-point (N k ).For each data point in this figure, we label the average atomic force for each atom along a certain direction (γ ∈ x, y, z) as F ave , which is computed by averaging 9 independent SDFT calculations with different sets of stochastic orbitals.In each SDFT calculation, we label the force acting on each atom along the γ direction as F sto .In this regard, the root-mean-square error (RMSE) can be evaluated via where N c = 9 is the number of independent SDFT runs, and N = 32 is the number of atoms with i being the index of atoms.
The number of stochastic orbitals is chosen from 128 to 11082 in Fig. 2(a), and the shifted k-point sampling is set to 2×2×2.Note that after symmetry analysis, the number of k points reduces to 4 after symmetry analysis.We fix the number of stochastic orbitals to be 512 in Fig. 2(b) and choose the shifted k-point samplings of , and (N sto × N k ) −1/2 .The numerical results are consistent with the discussion of statistical error in Section II A. Importantly, the results indicate that as more stochastic orbitals and a larger number of k-point sampling are employed, the stochastic errors can be systematically mitigated.

B. Compare SDFT and KSDFT Results
In order to select suitable numbers of k-points, KS orbitals, and stochastic orbitals for simulating warm dense matter under specific conditions, we use a 32-atom B system as an example, with a density of 2.46 g/cm 3 and a temperature of 17.23 eV. Figure 3 compares the forces on each B atom with different values for the mentioned parameters.First, we find that increasing the k-point sampling size from the Γ point to a shifted 2 × 2 × 2 kpoint sampling substantially reduces the RMSE, which is consistent across various values of N KS and N sto .The result is in line with the linear relationship shown in Fig. 2. Second, we note that the RMSE with N KS = 0 and N sto = 256 in Fig. 3(c) is 1.61 eV/ Å , while the RMSE with N KS = 128 and N sto = 128 shown in Fig. 3(e) is 0.596 eV/ Å .The former is considerably larger than the latter, suggesting that increasing the number of KS orbitals is more effective than using stochastic orbitals at a relatively low temperature (17.23 eV).Consequently, by choosing an adequate number of KS orbitals (N KS = 128) and stochastic orbitals (N sto = 128) along with the shifted 2 × 2 × 2 k-point sampling, we can achieve an RMSE as small as 0.302 eV/ Å .Furthermore, we examine the effects of these parameters on B systems at 86 and 350 eV, with the results shown in Figs.S1 and S2 of Supporting Information (SI), respectively.In conclusion, we find it reasonable to select the same parameters as in the B system at 17.23 eV.
Next, Table I presents a comparison of some key physical properties obtained from the SDFT and KSDFT methods.Specifically, we evaluate the total energy per atom (E), the pressure (P ), and the degree of ionization (α).∆ is the percentage difference between the results obtained from SDFT and the traditional KSDFT.We consider four systems, including two B systems at a temperature of 17.23 eV and densities of 2.46 and 12.3 g/cm 3 , as well as two C systems at a temperature of 21.54 eV and densities of 4.17 and 12.46 g/cm 3 .In both SDFT and KSDFT calculations, we choose the PBE 61 functional.Furthermore, a shifted 2 × 2 × 2 k-point sampling grid is adopted.
We study B systems with a cell containing 32 atoms.Additionally, we employ a norm-conserving pseudopotential for B with 3 valence electrons 65 , and the energy cutoff is 150 Ry.We respectively set the number of deterministic orbitals in KSDFT to be N KS = 992 and 400 for the B systems with density being 2.46 and 12.3 g/cm 3 .This setting ensures the occupation of electrons to be smaller than 10 −4 at the highest-energy orbital.On the other hand, we choose the number of deterministic orbitals to be N KS = 240 and the number of stochastic orbitals to be N sto = 120 in SDFT for the B systems regardless of their F a v e ( e V / Å ) online) Forces acting on each B atom as obtained from 9 independent SDFT calculations with different stochastic orbitals.The B system has a density of 2.46 g/cm 3 , and the temperature is 17.23 eV.For each calculation, the force acting on each atom along the γ direction is denoted as Fsto, and their average is Fave.NKS refers to the number of KS orbitals, and Nsto is the number of stochastic orbitals.Two sets of Monkhorst-Pack k-points are utilized, i.e., a 1 × 1 × 1 k-point mesh (the Γ k-point) and a shifted 2 × 2 × 2 k-point mesh.The root-mean-square error (RMSE) of forces between Fsto and Fave is computed via Eq. 13.The RMSE is obtained via the mentioned 9 independent SDFT calculations.
As shown in Table I, we have the following findings.First, it can be seen that the percentage difference in total energy (∆ of E) between SDFT and KSDFT is relatively small, being less than 0.02% for the B systems and 0.004% for the C systems.This indicates that SDFT provides a high-accuracy estimation of total energy when compared to the conventional KSDFT method.Second, the percentage difference in pressure (∆ of P ) between the two methods is smaller than 0.41% for B and 0.07% for C.This further supports the high accuracy of the SDFT method.Third, the ionization process of electrons plays a crucial role in determining the WDM equation of state 19,66,67 .This process can be represented by the degree of ionization, denoted as α.In practice, the Fermi energy of the system at 0 K is defined as µ.Consequently, the degree of ionization α at a finite temperature T can be defined as follows where N T,occ is the total number of occupied electrons below µ when the electrons follow the Fermi-Dirac distribution at tempearture T .In addition, N 0,occ is a special case of N T,occ when T = 0.The percentage difference in the degree of ionization (∆ of α) is found to be smaller than 0.009% for both B and C systems.In summary, all three properties, namely E, P , and α calculated by SDFT show excellent accuracy when compared to those from the traditional KSDFT.This demonstrates that SDFT is a reliable method for simulating hightemperature materials with first-principles accuracy.Fig. 4 further compares the forces acting on each atom of B (2.46 and 12.3 g/cm 3 ) and C (4.17 and 12.46 g/cm 3 ) obtained from both SDFT and KSDFT calculations.We find that the forces predicted by SDFT are in excellent agreement with those from KSDFT.For instance, the RMSE of forces is smaller than 0.05 eV/Å for both C systems.The largest RMSE occurs in the B system at 2.46 g/cm 3 , with a value of 0.103 eV/Å, which is relatively small compared to the magnitude of atomic forces (a few hundreds of eV/Å).Notably, we find the smallest RMSE is 0.004 eV/Å in the B system at 12.3 g/cm 3 .This is due to the fact that more electronic states of B are occupied by electrons at this condition, as demonstrated by the smaller degree of ionization of B (0.39) shown in Table .I.
Fig. 5 illustrates the DOS of B (2.46 and 12.3 g/cm 3 ) and C (4.17 and 12.46 g/cm 3 ).Besides the PBE 61 exchange-correlation functional, we also test the finitetemperature local density approximation functional, i.e., the corrKSDT functional as proposed by Karasiev et al. 64 .A Monkhorst-Pack 4 × 4 × 4 shifted k-point mesh is adopted in KSDFT calculations to yield the DOS of B and C.However, unlike the traditional KSDFT method, DOS in SDFT cannot be directly obtained from the eigenvalues of Ĥ.Instead, we evaluate the DOS from the SDFT method via the following formula Here, σ depicts the width of smearing.The DOS of SDFT with a shifted 2×2×2 k-point mesh converges for B with a density of 2.46 g/cm 3 when compared with KSDFT, although there are some deviations observed in the other three cases.Notably, the DOS predicted by SDFT using a 4 × 4 × 4 shifted k-point mesh shows excellent agreement with the KSDFT results for both B and C systems.By employing a shifted 4 × 4 × 4 k-point mesh, it is also observed that the DOS of corrKSDT 64 exhibits no significant differences when compared with the PBE 61 results, suggesting that the temperature effects in the XC functional have minimal impacts on our calculations.Overall, these findings indicate that the SDFT implemented in ABACUS is adequately accurate for simulating warm dense B and C systems.

C. High-Temperature Calculations by SDFT
Table II collects the three force components of each atom in the 32-atom B cell from 9 independent runs of SDFT with different stochastic orbitals.Four different (Color online) Density of states (DOS) for the B and C systems as calculated by the SDFT and KSDFT methods.The densities are selected as (a) 2.46 g/cm 3 and (b) 12.3 g/cm 3 for B, and (c) 4.17 g/cm 3 and (d) 12.64 g/cm 3 for C systems.The Fermi energy is set to 0. We use two sets of k-point sampling in KSDFT and SDFT calculations.In addition, we adopt two exchange-correlation functionals, including PBE 61 and corrKSDT 64 .TABLE II.Root-Mean-Square Error (RMSE) of forces acting on B atoms as obtained from 9 independent SDFT calculations with two sets of stochastic orbitals (128 and 256).The B system has a density of 2.46 g/cm 3 , and the temperatures are set to 17.23, 86, 350, and 1000 eV.NKS refers to the number of KS orbitals, and Nsto is the number of stochastic orbitals.Two sets of Monkhorst-Pack k-points are utilized, i.e., a 1 × 1 × 1 k-point mesh (include the Γ k-point) and a shifted 2 × 2 × 2 k-point mesh.The RMSE of forces is computed via Eq. 13

D. Radial Distribution Functions
Previous works have employed the traditional KSDFT coupling with BOMD to study WDM at relatively low  19 The SDFT calculations are performed with the PBE 61 and corrKSDT 64 exchange-correlation functionals.The number of B atoms is set to 32 in the first-principles molecular dynamics simulations.DPMD depicts the model trained by the traditional DP method 50 , while DPMD-T means the model trained by the TDDP method 51 .N is the number of B atoms in a cell, which ranges from 32 to 16384 in the deep-potential-based simulations.
temperatures.The examples include the shock Hugoniot curves 30 , the radial distribution function 11,15 , the ion-ion static structure factor 15 , and the ion-ion dynamic structure factor 15,55 , etc.However, most of these calculations are limited due to the high computational costs of dealing with electrons at extremely high temperatures.In order to substantially accelerate the calculations, we choose the DPMD method to perform BOMD calculations for warm sense B systems, and the training data are obtained from efficient SDFT calculations for warm dense B at temperatures of 86 and 350 eV.Note that we use the traditional DP method 50 to train the data at a temperature of 86 eV.However, we utilize the TDDP method 51 and include the electron temperature as an input parameter of the neural network to train the high-temperature data (350 eV), which shows a better performance than the traditional DP method at such a high temperature.
In detail, we perform SDFT-based BOMD simulations for a 32-atom B system with a density of 2.46 g/cm 3 .At the temperature of 86 eV, we adopt a Γ k-point mesh with the number of Kohn-Sham orbitals being N KS = 128 and the number of stochastic orbitals being N sto = 128.For a higher temperature of 350 eV, a 2×2×2 shifted k-point mesh is used with the number of stochastic orbitals being N sto = 256.We note that convergence with respect to the plane-wave energy cutoff and k-point mesh is examined to ensure the computational error of the total energy is within 0.01%.The BOMD simulations are performed in the NVT ensemble with the ion temperature controlled by the velocity-rescaling thermostat.The electrons and ions in the system are set to the same temperature.The time step is chosen according to ∆t , where T e is the temperature of electrons and ρ is the density.As a result, the time step is chosen to be 0.035 and 0.007 fs for simulations at 86 and 350 eV, respectively.In each BOMD trajectory, 4000 MD steps are performed.We then collect the atomic positions, the total energies E, the atomic forces F i of each atom i, as well as the virial tensors Ξ as the training data to generate DP models for B. Although stochastic DFT exhibits favorable scalability with increasing atoms 20 , previous research 51,[68][69][70] suggests that having 32 Boron atoms in the training dataset is enough to generate an accurate deep potential.The use of 32 Boron atoms to generate the training set is a choice that balances efficiency and accuracy.
In DPMD simulations, we adopt the NVT ensemble with the Nosé-Hoover thermostat 71,72 .We use the LAMMPS package 73 .The number of B atoms ranges from 32 to 16384.A time step of 0.07 fs is set for the system at 86 eV and 0.01 fs for the system at 350 eV.We perform 400000 steps of DPMD simulations to yield the radial distribution functions, the static structure factors, and the dynamic structure factors.Furthermore, one million MD steps are performed to compute the shear viscosity.
After the BOMD trajectories are generated, the radial distribution function can be evaluated according to where V is the cell volume, N is the number of atoms, r i and r j are atomic coordinates of atoms i and j, and ⟨• • • ⟩ means the the ensemble average.We plot g(r) of warm dense B with a density of 2.46 g/cm 3 at 86 and 350 eV in Fig. 6.We use Ext-FPMD 19 , SDFT with two different XC functionals (PBE and corrKSDT), and DPMD methods to perform BOMD simulations.We have the following findings.First, we find that the SDFT results are in excellent agreement with those obtained from Ext-FPMD.Second, there are no substantial differences between the PBE 61 and the finite-temperature XC functional corrKSDT 64 , which indicates that temperature effects in the exchange-correlation functional are not sig-nificant.Third, as expected, the g(r) is not smooth due to a limited number of MD steps (4000).However, by employing the DPMD method, we not only achieve accurate g(r) that agree well with the SDFT results but also obtain a smooth g(r), as a larger number of atoms (108 to 16384) and a longer trajectory (400000 steps) are considered.Importantly, size effects can be largely mitigated, as evidenced by the convergence of the g(r) at around 108 atoms.The ion-ion static structure factor S(q) measured from diffraction experiments 74,75 contains information regarding the spatial arrangement of particles in a material.

E. Ion-Ion Static Structure Factors
The formula of S(q) is where N is the number of atoms, while i, j denote atoms, and q is a wave vector.Here, we perform BOMD simulations on a 32-atom cell by SDFT with the PBE 61 and corrKSDT 64 exchange-correlation functionals.Moreover, we employ DPMD to calculate the S(q) for cells containing 32, 108, 256, 2048, and 16384 B atoms with a density of 2.46 g/cm 3 .The results for systems at 86 and 350 eV are illustrated in Figs.7(a) and (b), respectively.It is noteworthy that the data points of S(q) generated by SDFT exhibit oscillations due to the limited number of MD steps (4000).However, the DPMD simulations offer more converged results as they allow for a larger cell size with considerably more atoms and a substantially longer trajectory in BOMD simulations.Furthermore, with the use of larger cells in DPMD, we can obtain reasonable low-q information of S(q), which signifies the long-ranged order of systems.

F. Ion-Ion Dynamic Structure Factors
The collective dynamics of ionic density fluctuations can be characterized by the ion-ion dynamic structure factor S(q, ω), which is experimentally measurable 76 and plays a crucial role in investigating ion dynamics, including collective modes 77 , dissipation processes 78 , and others.In practice, S(q, ω) can be computed from the intermediate scattering function F (q, t) via Fourier transform with the formula being Here F (q, t) takes the form of where N is the number of ions, ρ(q, t) is defined as where r i (t) is the atomic coordinates for atom i at time t.
Figs. 8(a) and (b) illustrate the intermediate scattering function F (q, t) of warm dense B at 86 and 350 eV, respectively.Three wave vectors are chosen (q=0.51,1.02, and 2.50) while three sizes of cells are tested (256, 2048, and 16384 atoms).We find the 256-atom cell is large enough to converge F (q, t) for both temperatures, which is beyond the capabilities of SDFT-based BOMD simulations.
online) Intermediate scattering functions F (q, t) of B with a density of 2.46 g/cm 3 as calculated from DPMD trajectories.Three system sizes (256, 2048, and 16384 atoms) are adopted in DPMD simulations at two temperatures of (a) 86 eV and (b) 350 eV.Three wave vectors are chosen, i.e., q=0.51, 1.02, and 2.50 Å −1 .
Next, we obtain the ion-ion dynamics structure factors S(q, ω) of warm dense B by performing the Fourier transform of F (q, t).The results associated with wave vectors q being 0.51, 1.02, and 2.50 Å −1 are shown in Figs.9(a), (b), and (c), respectively.In each figure, two temperatures (86 and 350 eV) and three system sizes (256, 2048, and 16384 atoms) are adopted.For the wave vector, q=0.51 Å −1 , Fig. 9(a) shows the well-pronounced ionacoustic modes located at ω = 206.78meV for 86 eV and 486.80 meV for 350 eV.When q increases to 1.02 Å −1 in Fig. 9(b), the peak of S(q, ω) becomes lower, and its location turns to 324.95 meV for 86 eV and 616.53 meV for 350 eV.Notably, the ion-acoustic mode S(q, ω) disappears when q = 2.50 Å −1 , because the non-collective mode dominates at large q.The above results of S(q, ω) demonstrate that the DP method can predict the longranged structural and time correlation of WDM.For high temperatures up to hundreds of eV, there are experimental measurements of the ion-ion static structure factors and dynamic structure factors via X-ray Thomson scattering for materials, such as CH 79 and Be 80 .To the best of our knowledge, no experimental data are available for B at the temperatures of 86 and 350 eV.(Color online) Ion-ion dynamic structure factors S(q, ω) of B with a density of 2.46 g/cm 3 as computed from DPMD trajectories.Three system sizes (256, 2048, and 16384 atoms) are adopted.The wave vectors are chosen to be (a) q=0.51 Å −1 , (b) 1.02 Å −1 , and (c) 2.50 Å −1 .

G. Shear Viscosities
The shear viscosity is a crucial parameter in WDM studies, but obtaining a converged viscosity using traditional first-principles molecular dynamics is computationally expensive.However, this challenge can be significantly mitigated by employing the DPMD method with the training data from the SDFT method.One way to compute the shear viscosity η of WDM is using the Green-Kubo relations 81,82 , which links the shear viscosity to the integral of the stress auto-correlation function with the form of where V is the volume of the system, T is the temperature, k B is the Boltzmann constant, and P αβ (t) (αβ ∈ {xy, xz, yz}) is any of the three independent off-diagonal elements of the stress tensor at time t.The above formula can be used when DPMD trajectories are generated with the stress tensors.
The calculated stress auto-correlation functions of B at a density of 2.46 g/cm 3 and temperatures of 86 and 350 eV are displayed in Figs.10(a) and (b), respectively.In practice, the computed shear viscosity may be affected by the system size and the trajectory length of molecular dynamics simulations.Therefore, we choose seven different system sizes with the number of atoms per cell ranging from 32 to 16384.During DPMD simulations, each system is first relaxed for 50000 steps.Next, 1 million steps of MD simulations are performed to calculate the stress auto-correlation function.In detail, the trajectory length is 70 ps for 86 eV and 10 ps for 350 eV.We take values from 0.105 to 0.305 ps (0.05 to 0.1 ps) to compute the averaged shear viscosity for the system at 86 eV (350 eV), and the predicted values are shown in Fig. 10(c) with small error bars.As a result, the obtained shear viscosity of B varies from 10.3 to 12.3 mPa • s at 86 eV and from 44.4 to 47.3 mPa • s at 350 eV.The above results show no substantial finite-size effects for the shear viscosity of B, which is consistent with previous conclusions [83][84][85] .
There are other formulas that can predict the shear viscosity of materials.For example, we notice that a recent work proposes an extended random-walk shieldingpotential viscosity model (ext-RWSP-VM) 86,87 to elevate the shear viscosity of materials in WDM and HDP states.The viscosity is evaluated by the formula of where d is the collision diameter introduced by hardsphere concept, and I is a quantity that is relevant to T .According to the ext-RWSP-VM method, we obtain the viscosities of B to be 12.8 and 47.8 mPa • s for temperatures of 86 and 350 eV, respectively.More details can be found in SI.Interestingly, the data are close to our first-principles results.In addition, we find the shear viscosity of plasma can also be described by the approximate formula 88 of where m is the atomic mass and σ o is the total collision cross section (∼ 10 −20 m 2 ).Thus, the estimated shear viscosities of B are around ∼ 18.7 and 37.8 mPa • s for temperatures of 86 and 350 eV, respectively.It should be noted that σ o in the approximate formula is assumed to be a constant; however, it is related to the relative velocity between atoms 88 .As the relative velocity increases, the interaction time decreases, leading to a reduced probability of collisions occurring.In other words, σ o decreases with increasing temperature.We find that the calculated shear viscosities from DPMD are also consistent with the approximated values obtained using Eq.23.

IV. CONCLUSIONS
Simulating WDM with first-principles accuracy had long been challenging due to the existence of the partial occupation of a large number of high-energy KS eigenstates and the resulting limitations in the time and space scales.Our work suggested that the advent of the SDFT method and machine-learningbased molecular dynamics can be of great help in overcoming the difficulties.The SDFT method described in this work had been implemented with the plane-wave basis set and the k-point sampling method, which was enabled in the ABACUS package (https://github.com/deepmodeling/abacus-develop).In this work, we validated the SDFT-based BOMD method by performing a series of tests for warm dense B and C.
By combining SDFT with the DP method, we substantially extended the time and space scales of simulating warm dense B and reduced the finite size effect.Besides, we studied the structural properties, dynamic properties, and transport coefficients, such as radial distribution functions, static structure factors, ion-ion dynamic structure factors, and shear viscosities.This work validated combining stochastic density functional theory with machine learning techniques to study high-temperature systems.We also offered new insights into the properties of warm dense matter.In future work, we intend to explore the generation of training data with a larger number of atoms.Future research may further refine these methods and expand their applicability to other materials and temperature ranges.
distance equals the Debye length where ε 0 is vacuum permittivity, n e is the number density of the electrons, and Here, we assume that the ionization is uniform, so z * = Z = 4.30788 for the temperature of 86 eV and z * = 4.93598 for 350 eV.The average ionization Z = Z × α, where α is the degree of ionization calculated by Eq. ( 14) in the main article, and Z is the nuclear charge.
The collision diameter takes the formula of where R = r 0 − a, and b m = r 2 0 − 2r 0 a.

FIG. 1 .
FIG. 1. (Color online) Workflow of simulating WDM with the SDFT and DPMD methods.(a) We employ stochastic orbitals in SDFT to perform MD simulations on smaller systems (32 atoms in a bulk B) and collect initial training data, which include atomic positions, energies, forces, and virial tensors.(b) We utilize the gathered training data to construct a DP model with the temperature-dependent DPMD model.The deep neural network contains both embedding and fitting networks.(c) The new model enables us to perform simulations on large systems (16384 atoms) and at extremely high temperatures (350 eV).(d) Several physical quantities such as radial distribution functions, static structure factors, dynamic structure factors, and shear viscosities can be calculated, and the data are converged with large systems and long trajectories.

FIG. 2 .
FIG. 2. (Color online) Root-mean-square error (RMSE) of atomic forces arises from the SDFT calculations for B. The temperature is set to 350 eV and the density is set to 2.46 g/cm 3 .The RMSE is evaluated with respect to (a) the number of stochastic orbitals Nsto, (b) the number of k points in the Brillouin zone N k , and (c) the product of Nsto and N k .For each data point, the RMSE is obtained via 9 independent SDFT calculations with different sets of stochastic orbitals.
and 5 × 5 × 5 (65); here the number in the parentheses denotes the number of kpoints after symmetry analysis.As expected, Figs.2(a), (b), and (c) respectively show that the RMSE of forces acting on atoms exhibits linear behavior with respect to N

FIG. 4 .
FIG. 4. (Color online) (a) Comparison of forces acting on each B atom in a 32-atom cell with densities being 2.46 and 12.3 g/cm 3 at a temperature of 17.23 eV.(b) Comparisonof forces acting on each C atom with densities being 4.17 and 12.46 g/cm 3 at a temperature of 21.54 eV.In the SDFT (KSDFT) calculation, we label the force acting on each atom along the γ direction (γ ∈ x, y, z) as Fsto (FKS).The rootmean-square error (RMSE) of forces between Fsto and FKS is listed.
FIG. 5.(Color online) Density of states (DOS) for the B and C systems as calculated by the SDFT and KSDFT methods.The densities are selected as (a) 2.46 g/cm 3 and (b) 12.3 g/cm 3 for B, and (c) 4.17 g/cm 3 and (d) 12.64 g/cm 3 for C systems.The Fermi energy is set to 0. We use two sets of k-point sampling in KSDFT and SDFT calculations.In addition, we adopt two exchange-correlation functionals, including PBE 61 and corrKSDT64 .

FIG. 6 .
FIG.6.(Color online) Radial distribution functions g(r) of B systems with a density of 2.46 g/cm 3 at temperatures of (a) 86 eV and (b) 350 eV.The g(r) obtained by Ext-FPMD with the LDA functional comes from Blanchet et al.19The SDFT calculations are performed with the PBE61 and corrKSDT64 exchange-correlation functionals.The number of B atoms is set to 32 in the first-principles molecular dynamics simulations.DPMD depicts the model trained by the traditional DP method50  , while DPMD-T means the model trained by the TDDP method51  .N is the number of B atoms in a cell, which ranges from 32 to 16384 in the deep-potential-based simulations.
FIG. 7.(Color online) Static structure factors S(q) of B with a density of 2.46 g/cm 3 at temperatures of (a) 86 eV and (b) 350 eV.The SDFT calculations are performed with the PBE 61 and corrKSDT 64 exchange-correlation functionals.The number of B atoms is set to 32 in the first-principles molecular dynamics simulations.DPMD depicts the model trained by the traditional DP method 50 , while DPMD-T indicates the model trained by the TDDP method51  .N is the number of B atoms in a cell ranging from 32 to 16384 in the deep-potential-based simulations.
FIG. S1. (Color online) Forces on each atom of B with a density of 2.46 g/cm 3 for 86 eV calculated by SDFT.F sto is the force calculated by SDFT with nine different random seeds to generate stochastic orbitals, and their average is F ave .N KS is the number of KS orbitals and N sto is the number of stochastic orbitals.1 × 1 × 1 k-point means the Γ k-point and 2 × 2 × 2 k-point means × 2 × 2 shifted k-point.RMSE is the root mean square error of forces between F sto and F ave .
FIG.S3.(Color online) Forces acting on each B atom as obtained from 9 independent SDFT calculations with different stochastic orbitals.The B system has a density of 2.46 g/cm 3 , and the temperatures are set to 17.23 eV, 86 eV, 350 eV, and 1000 eV.For each calculation, the force acting on each atom along the γ direction is denoted as F sto , and their average is F ave .N KS refers to the number of KS orbitals, and N sto is the number of stochastic orbitals.Two sets of Monkhorst-Pack k-points are utilized, i.e., a 1 × 1 × 1 k-point mesh (the Γ k-point) and a shifted 2 × 2 × 2 k-point mesh.RMSE is the root mean square error of forces between F sto and F ave .
FIG. S4. (Color online) Changes of (a) temperature (in K), (b) total energy (in Ha), and (c) pressure (in kbar) with respect to MD steps when simulating B systems at temperatures of 86 and 350 eV with the PBE ? and corrKSDT ?exchange-correlation functionals.

TABLE I .
Comparison of the total energy per atom (E in eV), the pressure (P in GPa), and the degree of ionization (α) for B and C systems as obtained from the SDFT and the traditional KSDFT methods.Four systems are chosen, i.e., two B systems at a temperature of 17.23 eV with densities being 2.46 and 12.3 g/cm 3 , and two C systems with densities being 4.17 and 12.46 g/cm 3 at a temperature of 21.54 eV.∆ denotes the percentage difference between the results obtained by SDT and KSDFT. .
temperatures are chosen, i.e.,17.23,86,350, and 1000 eV.Furthermore, we select 1 × 1 × 1 and 2 × 2 × 2 shifted k-point samplings for each temperature and evaluate the corresponding RMSE.At each condition, a set of average forces F ave are calculated according to Eq. 13.One can refer to Fig.S3of SI for more details.For temperatures of 17.23 and 86 eV, we set N KS = 128 and N sto = 128;