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

Tao Chen; Qianrui Liu; Yu Liu; Liang Sun; Mohan Chen

doi:doi:10.1063/5.0163303

Matter and Radiation at Extremes, 2024, 9 (1): 015604, Published Online: Mar. 27, 2024

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

论文大纲

Tao Chen Qianrui Liu Yu Liu Liang Sun Mohan Chen ^a)

Author Affiliations

HEDPS, CAPT, College of Engineering and School of Physics, Peking University, Beijing 100871, People’s Republic of China

Abstract

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 this limitation. Recently, SDFT and the related mixed stochastic–deterministic density functional theory, based on a plane-wave basis set, have been implemented in the first-principles electronic structure software ABACUS [Q. Liu and M. Chen, Phys. Rev. B 106, 125132 (2022)]. In this study, we combine SDFT with the Born–Oppenheimer molecular dynamics 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. Subsequently, we compute and analyze the structural properties, dynamic properties, and transport coefficients of warm dense matter.

1 I. INTRODUCTION

Understanding the behavior of materials under extremely high-temperature conditions, such as warm dense matter (WDM) and hot dense plasma (HDP), is essential in a number of physical and astrophysical contexts, such as the interiors of giant planets,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 experimental5 and theoretical4 perspectives, it has been recognized that an adequate quantum-mechanical description of electrons is essential for theoretical calculations to have the necessary predictive power. In fact, various first-principles approaches based on different levels of approximations are available to address this issue. The methods most commonly employed for this task are Kohn–Sham density functional theory (KSDFT),6,7 path-integral Monte Carlo (PIMC),8–11 orbital-free density functional theory (OFDFT),12–15 extended first-principles molecular dynamics (Ext-FPMD),16–19 and stochastic density functional theory (SDFT).20–24

KSDFT with the Mermin finite-temperature density-functional approach25 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, since the number of electronic states that need to be considered is proportional to T^3/2, rendering the KSDFT method unfeasible.26–29 By contrast, the cost of PIMC calculations decreases significantly at higher temperatures.8–11 In practice, although combining KSDFT with PIMC has been successfully applied to study the equations of state for low-Z elements,30 PIMC is severely limited at lower temperatures.31 Compared with KSDFT, the OFDFT method avoids the need to diagonalize the wave functions and is more efficient.12–15 However, the applications of OFDFT are still limited by insufficient accuracy in the kinetic energy density functional.32

The Ext-FPMD method has been proposed to evaluate systems at extremely high temperatures with first-principles accuracy.16–19 This method treats the wave functions of high-energy electrons analytically as plane waves, 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, owing to a 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 limitation imposed by partial occupation of a large number of high-energy KS eigenstates. Moreover, the statistical errors can be effectively reduced when a greater number of stochastic orbitals are included or larger systems are utilized. In addition, 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 a plane-wave basis set, can be used with the k-point sampling method and have recently been implemented in the first-principles package ABACUS (https://github.com/deepmodeling/abacus-develop).34–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 have attracted considerable attention.37–48 In particular, deep neural networks are often adopted to learn the potential energy surfaces that are generated 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, since they effectively bypass the need to solve quantum-mechanics-based equations. Among them, the deep potential molecular dynamics (DPMD) method49–52 achieves high performance53,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 properties15,55 and ion–ion dynamical structure factor55 of warm dense aluminum, thermal transport by electrons and ions in warm dense aluminum,56 and the principal Hugoniot curves of warm dense beryllium.51

To summarize, two main challenges exist for large-scale first-principles simulations of WDM. First, the partial occupation of a large number of high-energy 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 molecular dynamics trajectory is obtained.

In this work, we first validate the accuracy of SDFT by analyzing the statistical errors from the stochastic orbitals and compare our results with those from the traditional KSDFT 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³ with two different temperatures (86 and 350 eV); the training data are obtained from SDFT-based 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. This workflow is shown in Fig. 1.

Workflow of the simulation of WDM with the SDFT and DPMD methods. (a) Stochastic orbitals are employed in SDFT to perform molecular dynamics simulations on smaller systems (32 atoms in a bulk B) and initial training data are collected that include atomic positions, energies, forces, and virial tensors. (b) The gathered training data are used 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 simulations to be performed on large systems (16 384 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. 1. Workflow of the simulation of WDM with the SDFT and DPMD methods. (a) Stochastic orbitals are employed in SDFT to perform molecular dynamics simulations on smaller systems (32 atoms in a bulk B) and initial training data are collected that include atomic positions, energies, forces, and virial tensors. (b) The gathered training data are used 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 simulations to be performed on large systems (16 384 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.

下载图片查看所有图片

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. In Sec. II, we describe the computational methods utilized in this work. In Sec. III, we discuss the results obtained from SDFT and the DP model. Finally, we summarize our results in Sec. IV.

2 II. COMPUTATIONAL METHODS

2.1 A. Stochastic and mixed stochastic–deterministic density functional theory

In the KSDFT framework,6,7 the single-particle DFT Hamiltonian is defined as ${\hat{H}}_{DFT} = - \nabla^{2} / 2 + {\hat{V}}_{Hartree} + {\hat{V}}_{xc} + {\hat{V}}_{ext},$ (1)where the first term is the kinetic operator, ${\hat{V}}_{Hartree}$ is the Hartree potential, ${\hat{V}}_{xc}$ is the exchange-correlation (XC) potential, and ${\hat{V}}_{ext}$ is the potential of interactions between the electrons and the nuclei as well as other external fields. Within finite-temperature Mermin–Kohn–Sham theory,25 the electron density is given by $\begin{align} n (r) & = 2 Tr [f ({\hat{H}}_{DFT}, μ, T) δ (\hat{r} - r)] \\ = 2 \sum_{k = 1}^{N_{k}} ω_{k} \sum_{i = 1}^{N_{KS}} f (ε_{i k}, μ, T) | ϕ_{i k} (r) |^{2}, \end{align}$ (2)where $\hat{r}$ is the position operator of the electron and 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 $f (ε_{i k}, μ, T) = 1 / [e^{(ε_{i k} - μ) / T} + 1]$ , where μ is the chemical potential. Here, ϕ_ik(r) and ɛ_ik are the eigenfunction and eigenvalue, respectively, of the self-consistent KS Hamiltonian ${\hat{H}}_{DFT} ϕ_{i k} (r) = ε_{i k} ϕ_{i k} (r) .$ (3)In practice, solving ϕ_ik(r) and ɛ_ik is costly at high temperatures, since the required 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 SDFT20,21 can be defined as $⟨ ψ_{j} ∣ χ_{n k} ⟩ = \frac{1}{\sqrt{N_{sto}}} \exp (i 2 π θ_{j}^{n k}),$ (4)which satisfies $\lim_{N_{sto} \to + \infty} \sum_{n = 1}^{N_{sto}} | χ_{n k} 〉 〈 χ_{n k} | = \hat{I},$ (5)where $θ_{j}^{n k}$ 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 ${\tilde{χ}}_{n k}$ are used: $| {\tilde{χ}}_{n k} 〉 = | χ_{n k} 〉 - \sum_{i = 1}^{N_{KS}} ⟨ ϕ_{i k} ∣ χ_{n k} ⟩ | ϕ_{i k} 〉 .$ (6)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 $\lim_{N_{sto} \to + \infty} \sum_{n = 1}^{N_{sto}} | {\tilde{χ}}_{n k} 〉 〈 {\tilde{χ}}_{n k} | + \sum_{i = 1}^{N_{KS}} | ϕ_{i k} 〉 〈 ϕ_{i k} | = \hat{I} .$ (7)The electron density is then given by $\begin{align} n (r) & = 2 \sum_{k = 1}^{N_{k}} ω_{k} [\sum_{n = 1}^{N_{sto}} {|⟨r |\sqrt{f ({\hat{H}}_{DFT}, μ, T)}| {\tilde{χ}}_{n k}⟩|}^{2} \\ + \sum_{i = 1}^{N_{KS}} f (ε_{i k}, μ, T) | ϕ_{i k} (r) |^{2}], \end{align}$ (8)where $\sqrt{f ({\hat{H}}_{DFT}, μ, T)}$ is calculated by the Chebyshev expansion.57 If N_sto = 0 or N_KS = 0, the MDFT method changes to the standard KSDFT or SDFT method, 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,22,23,36 the error caused by the stochastic orbitals is proportional to $1 / \sqrt{N_{sto}}$ . However, when periodic boundary conditions (PBCs) with k-point sampling are considered and each k-point has N_sto stochastic orbitals, the resulting σ_s is proportional to $1 / \sqrt{N_{k} N_{sto}}$ , suggesting that the use of a greater number of k-points can reduce the stochastic errors. To evaluate the accuracy of SDFT, we perform both SDFT and KSDFT calculations for B and C bulk systems under extreme conditions.

In the first-principles molecular dynamics 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 $F_{i} = - \frac{\partial E}{\partial R_{i}},$ (9)where E is the sum of the electrons’ energy and the ion–ion repulsion energy, and R_i is the position of atom i. By utilizing a 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: $F_{i} = F_{i}^{L} + F_{i}^{NL} + F_{i}^{II},$ (10)where $F_{i}^{L}$ is the local pseudopotential force term, $F_{i}^{NL}$ is the nonlocal pseudopotential force term, and $F_{i}^{II}$ is the Ewald force term originating from ion–ion interactions. Furthermore, the stress is defined as $\begin{align} σ_{α β} & = - \frac{1}{V} \frac{\partial E}{\partial ϵ_{α β}} \\ = σ_{α β}^{T} + σ_{α β}^{L} + σ_{α β}^{NL} + σ_{α β}^{Hartree} + σ_{α β}^{xc} + σ_{α β}^{II}, \end{align}$ (11)where ϵ_αβ is the strain, with spatial coordinates α and β. The kinetic energy term of the electrons is $σ_{α β}^{T}$ . $σ_{α β}^{L}$ and $σ_{α β}^{NL}$ are the local and nonlocal pseudopotential terms, respectively. $σ_{α β}^{Hartree}$ is the Hartree term, and $σ_{α β}^{xc}$ is the XC term. The Ewald term is $σ_{α β}^{II}$ . Further details of the implementation of the total energy, total free energy, forces, and stresses within the framework of SDFT in ABACUS can be found in Ref. 36. Note that the SDFT method employed in this work still uses a plane wave basis, and so it is still computationally demanding for a system consisting of a few hundred atoms.

2.2 B. Deep potential molecular dynamics

In the DPMD method,49,50 the total energy E of a system is expressed as a sum of atomic contributions, i.e., E = ∑_iE_i, where the energy E_i from atom i depends on an environment matrix $R_{i}$ , which includes the information of those atoms neighboring atom i within a cutoff radius. The DP model maps $R_{i}$ via an embedding neural network to a symmetry-preserving descriptor and then maps this descriptor 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 $L (p_{ϵ}, p_{f}, p_{ξ}) = p_{ϵ} Δ ϵ^{2} + \frac{p_{f}}{3 N} \sum_{i} | Δ F_{i} |^{2} + \frac{p_{ξ}}{9} ‖ Δ ξ ‖^{2},$ (12)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 Δ is the difference between the training data and the results predicted by the DP model. Here, p_ϵ, p_f, and p_ξ are tunable prefactors. The stochastic gradient descent scheme Adam58 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)59 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) method51 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 numbers 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 Å 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 with both DP models to accelerate the DPMD simulations, as described in the literature.60

3 III. RESULTS AND DISCUSSION

3.1 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³ at a temperature of 350 eV. In addition, we employ the PBE59 XC functional.

We note that at temperatures of 86 eV and higher, the pseudopotential of B is generated by the ONCVPSP61 method with all of its five 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 energy cutoffs of 180, 240, and 300 Ry for temperatures of 86, 350, and 1000 eV, respectively.

Figure 2 shows the root-mean-square error (RMSE) of SDFT vs the number of stochastic orbitals N_sto and the number of k-points N_k.

Root-mean-square error (RMSE) of atomic forces arising from the SDFT calculations for B. The temperature is set to 350 eV and the density to 2.46 g/cm3. The RMSE is evaluated with respect to (a) the number of stochastic orbitals Nsto, (b) the number of k-points in the Brillouin zone Nk, and (c) the product of Nsto and Nk. For each data point, the RMSE is obtained via nine independent SDFT calculations with different sets of stochastic orbitals.

Fig. 2. Root-mean-square error (RMSE) of atomic forces arising from the SDFT calculations for B. The temperature is set to 350 eV and the density to 2.46 g/cm³. The RMSE is evaluated with respect to (a) the number of stochastic orbitals N_sto, (b) the number of k-points in the Brillouin zone N_k, and (c) the product of N_sto and N_k. For each data point, the RMSE is obtained via nine independent SDFT calculations with different sets of stochastic orbitals.

下载图片查看所有图片

For each data point in this figure, we denote the average atomic force for each atom along a certain direction γ ∈ x, y, z by F_ave, which is computed by averaging nine independent SDFT calculations with different sets of stochastic orbitals. In each SDFT calculation, we denote the force acting on each atom along the γ direction by F_sto. In this case, the RMSE can be evaluated as $R M S E = \sqrt{\frac{1}{3 N_{c} N} \sum_{c = 1}^{N_{c}} \sum_{i = 1}^{N} \sum_{γ} {(F_{sto} - F_{ave})}^{2}},$ (13)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 11 082 in Fig. 2(a), and the shifted k-point sampling is set to 2 × 2 × 2. Note that the number of k-points is reduced to four after symmetry analysis.

We fix the number of stochastic orbitals to be 512 in Fig. 2(b) and choose shifted k-point samplings of 1 × 1 × 1 (1), 2 × 2 × 2 (4), 3 × 3 × 3 (14), 4 × 4 × 4 (32), and 5 × 5 × 5 (63); here the numbers in parentheses are the numbers of k-points after symmetry analysis.

As expected, Figs. 2(a)–2(c) respectively show that the RMSE of forces acting on atoms exhibits linear behavior with respect to $N_{sto}^{- 1 / 2}$ , $N_{k}^{- 1 / 2}$ , and ${(N_{sto} N_{k})}^{- 1 / 2}$ . The numerical results are consistent with the discussion of statistical error in Sec. II A. Importantly, the results indicate that as more stochastic orbitals and a larger number of k-point samplings are employed, the stochastic errors can be systematically mitigated.

3.2 B. Comparison of SDFT and KSDFT results

To select suitable numbers of k-points, KS orbitals, and stochastic orbitals for simulating WDM under specific conditions, we use a 32-atom B system as an example, with a density of 2.46 g/cm³ 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 k-point 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.

Forces acting on each B atom as obtained from nine independent SDFT calculations with different stochastic orbitals. The B system has a density of 2.46 g/cm3, and the temperature is 17.23 eV. For each calculation, the force acting on each atom along the γ direction is denoted by Fsto, and their average is Fave. NKS and Nsto are the numbers of KS and stochastic orbitals, respectively. Two sets of Monkhorst–Pack k-points are utilized: a 1 × 1 × 1 k-point mesh (the Γ k-point) and a shifted 2 × 2 × 2 k-point mesh. The RMSE of forces between Fsto and Fave is computed from Eq. (13). The RMSE is obtained via the mentioned nine independent SDFT calculations.

Fig. 3. Forces acting on each B atom as obtained from nine independent SDFT calculations with different stochastic orbitals. The B system has a density of 2.46 g/cm³, and the temperature is 17.23 eV. For each calculation, the force acting on each atom along the γ direction is denoted by F_sto, and their average is F_ave. N_KS and N_sto are the numbers of KS and stochastic orbitals, respectively. Two sets of Monkhorst–Pack k-points are utilized: a 1 × 1 × 1 k-point mesh (the Γ k-point) and a shifted 2 × 2 × 2 k-point mesh. The RMSE of forces between F_sto and F_ave is computed from Eq. (13). The RMSE is obtained via the mentioned nine independent SDFT calculations.

下载图片查看所有图片

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 (supplementary material), respectively. In conclusion, we find it reasonable to select the same parameters as in the B system at 17.23 eV.

Table I presents a comparison of some key physical properties obtained from the SDFT and KSDFT methods, namely, the total energy per atom E, the pressure P, and the degree of ionization α. For each property, we also show the percentage difference Δ between the results obtained from SDFT and the traditional KSDFT. We consider four systems, namely, two B systems at a temperature of 17.23 eV and densities of 2.46 and 12.3 g/cm³, and two C systems at a temperature of 21.54 eV and densities of 4.17 and 12.46 g/cm³. In both the SDFT and KSDFT calculations, we choose the PBE functional.59 Furthermore, we adopt a shifted 2 × 2 × 2 k-point sampling grid.

Table 1. Comparison of the total energy per atom E, the pressure P, and the degree of ionization α for B and C systems as obtained from the SDFT and traditional KSDFT methods. Four systems are chosen: two B systems at a temperature of 17.23 eV with densities of 2.46 and 12.3 g/cm³, and two C systems with densities of 4.17 and 12.46 g/cm³ at a temperature of 21.54 eV. Δ denotes the percentage difference between the results obtained by SDFT and KSDFT.

		E (eV)	P (GPa)	α
B	KSDFT	−153.683 171	849.112	0.476 617
2.46 g/cm³	SDFT	−153.707 884	852.524	0.476 659
17.23 eV	Δ	0.0161%	0.4019%	0.0088%
B	KSDFT	−66.923 177	8730.010	0.390 341
12.3 g/cm³	SDFT	−66.924 235	8730.829	0.390 344
17.23 eV	Δ	0.0016%	0.0094%	0.0008%
C	KSDFT	−263.586 543	2018.380	0.505 693
4.17 g/cm³	SDFT	−263.595 810	2019.780	0.505 707
21.54 eV	Δ	0.0035%	0.0693%	0.0028%
C	KSDFT	−190.692 843	9168.050	0.440 747
12.46 g/cm³	SDFT	−190.698 535	9171.552	0.440 757
21.54 eV	Δ	0.0030%	0.0382%	0.0023%

查看所有表

We study B systems with a cell containing 32 atoms. Additionally, we employ a norm-conserving pseudopotential for B with three valence electrons62 and with an energy cutoff of 150 Ry. We set the numbers of deterministic orbitals in KSDFT to be N_KS = 992 and 400 for the B systems with densities of 2.46 and 12.3 g/cm³, respectively. These settings ensure that the occupation of electrons is smaller than 10⁻⁴ 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 densities. For the C systems, a norm-conserving pseudopotential with four valence electrons is employed,62 and the energy cutoff is set to 160 Ry. In the KSDFT calculations, we take 8 atoms per cell and N_KS = 350 for a density of 4.17 g/cm³ and 32 atoms per cell and N_KS = 520 for a density of 12.46 g/cm³. In the SDFT calculations, we adopt N_KS = 120 and 240 for densities of 4.17 and 12.46 g/cm³, respectively, with N_sto = 120 for both cases.

The results in Table I reveal the following.

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 with 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,63,64 This process can be represented by the degree of ionization α. 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: $α = 1 - \frac{N_{T, o c c}}{N_{0, o c c}},$ (14)where N_T,occ is the total number of occupied electrons below μ when the electrons follow the Fermi–Dirac distribution at temperature T. N_0,occ is the 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 E, P, and α calculated by SDFT show excellent accuracy when compared with those from traditional KSDFT. This demonstrates that SDFT is a reliable method for simulating high-temperature materials with first-principles accuracy.

Figure 4 further compares the forces acting on each atom of B (2.46 and 12.3 g/cm³) and C (4.17 and 12.46 g/cm³) obtained from both SDFT and KSDFT calculations.

(a) Comparison of forces acting on each B atom in a 32-atom cell for densities of 2.46 and 12.3 g/cm3 at a temperature of 17.23 eV. (b) Comparison of forces acting on each C atom for densities of 4.17 and 12.46 g/cm3 at a temperature of 21.54 eV. In the SDFT and KSDFT calculations, we denote the forces acting on each atom along the γ direction (γ ∈ x, y, z) by Fsto and FKS, respectively. The RMSE of forces between Fsto and FKS is also shown.

Fig. 4. (a) Comparison of forces acting on each B atom in a 32-atom cell for densities of 2.46 and 12.3 g/cm³ at a temperature of 17.23 eV. (b) Comparison of forces acting on each C atom for densities of 4.17 and 12.46 g/cm³ at a temperature of 21.54 eV. In the SDFT and KSDFT calculations, we denote the forces acting on each atom along the γ direction (γ ∈ x, y, z) by F_sto and F_KS, respectively. The RMSE of forces between F_sto and F_KS is also shown.

下载图片查看所有图片

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³, with a value of 0.103 eV/Å, which is relatively small compared with 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³. This is due to the fact that more electronic states of B are occupied by electrons under this condition, as demonstrated by the smaller degree of ionization of B (0.39) shown in Table I.

Figure 5 illustrates the density of states (DOS) of B (2.46 and 12.3 g/cm³) and C (4.17 and 12.46 g/cm³). Besides the PBE59 XC functional, we also test the finite-temperature local density approximation functional, i.e., the corrKSDT functional as proposed by Karasiev et al.65

DOS for the B and C systems as calculated by the SDFT and KSDFT methods. The densities are selected as (a) 2.46 g/cm3 and (b) 12.3 g/cm3 for the B system, and (c) 4.17 g/cm3 and (d) 12.64 g/cm3 for the C system. The Fermi energy is set to zero. We use two sets of k-point sampling in the KSDFT and SDFT calculations. In addition, we adopt two XC functionals, namely, PBE59 and corrKSDT.65

Fig. 5. 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³ and (b) 12.3 g/cm³ for the B system, and (c) 4.17 g/cm³ and (d) 12.64 g/cm³ for the C system. The Fermi energy is set to zero. We use two sets of k-point sampling in the KSDFT and SDFT calculations. In addition, we adopt two XC functionals, namely, PBE⁵⁹ and corrKSDT.⁶⁵

下载图片查看所有图片

A Monkhorst–Pack 4 × 4 × 4 shifted k-point mesh is adopted in the KSDFT calculations to yield the DOS of B and C. However, unlike the traditional KSDFT method, the DOS in SDFT cannot be directly obtained from the eigenvalues of $\hat{H}$ . Instead, we evaluate the DOS from the SDFT method via the following formula: $g (E) = 2 Tr \{\frac{1}{\sqrt{2 π} σ} \exp [- \frac{{(E - \hat{H})}^{2}}{2 σ^{2}}]\} .$ (15)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³ 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 the B and C systems.

By employing a shifted 4 × 4 × 4 k-point mesh, it is also observed that the DOS of corrKSDT65 exhibits no significant differences when compared with the PBE59 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.

3.3 C. High-temperature calculations by SDFT

Table II collects the three force components of each atom in the 32-atom B cell from nine independent runs of SDFT with different stochastic orbitals. Four different temperatures are chosen, namely, 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. Under each condition, a set of average forces F_ave are calculated according to Eq. (13). Further details can be found in Fig. S3 (supplementary material). For temperatures of 17.23 and 86 eV, we set N_KS = 128 and N_sto = 128; at higher temperatures of 350 and 1000 eV, we find it more effective to use stochastic orbitals than the KS orbitals. As a result, we do not choose the Kohn–Sham orbitals (N_KS = 0) but set all orbitals to be stochastic orbitals (N_sto = 256). According to our tests, the RMSE of the atomic force smaller than 3.3 eV/Å is enough for FPMD simulations of WDM B. Therefore, we employ the Γ k-point with 128 KS orbitals and 128 stochastic orbitals for 86 eV and the 2 × 2 × 2 shifted k-point with 256 stochastic orbitals for 350 eV to perform FPMD.

Table 2. RMSE of forces acting on B atoms as obtained from nine independent SDFT calculations with two sets of stochastic orbitals (128 and 256). The B system has a density of 2.46 g/cm³, and the temperatures are set to 17.23, 86, 350, and 1000 eV. N_KS is the number of KS orbitals, and N_sto is the number of stochastic orbitals. Two sets of Monkhorst–Pack k-points are utilized: 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 from Eq. (13).

T (eV)	N_KS	N_sto	k-points	RMSE (eV/Å)
17.23	128	128	1 × 1 × 1	0.596
	128	128	2 × 2 × 2	0.302
86	128	128	1 × 1 × 1	3.26
	128	128	2 × 2 × 2	1.54
350	0	256	1 × 1 × 1	6.40
	0	256	2 × 2 × 2	3.26
1000	0	256	1 × 1 × 1	2.68
	0	256	2 × 2 × 2	1.30

查看所有表

3.4 D. Radial distribution functions

Previous works have employed the traditional KSDFT coupling with BOMD to study WDM at relatively low temperatures. Examples include 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 However, most of these calculations have been limited by the high computational costs incurred when dealing with electrons at extremely high temperatures.

To substantially accelerate the calculations, we choose the DPMD method to perform BOMD calculations for warm dense 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 method50 to train the data at a temperature of 86 eV. However, we utilize the TDDP method51 and include the electron temperature as an input parameter of the neural network to train the high-temperature data (350 eV), since this method exhibits better performance than the traditional DP method at such a high temperature.

Specifically, we perform SDFT-based BOMD simulations for a 32-atom B system with a density of 2.46 g/cm³. At the temperature of 86 eV, we adopt a Γ k-point mesh with N_KS = 128 Kohn–Sham orbitals and N_sto = 128 stochastic orbitals.

At the higher temperature of 350 eV, a 2 × 2 × 2 shifted k-point mesh is used, with N_sto = 256 stochastic orbitals. Note that convergence with respect to the plane-wave energy cutoff and k-point mesh is examined to ensure tht 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 \sim 1 / (ρ^{1 / 3} T_{e}^{1 / 2})$ , where T_e is the temperature of the electrons and ρ is the density. Consequently, 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 molecular dynamics 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 number of atoms,20 previous research51,66–68 suggests that having 32 B atoms in the training dataset is enough to generate an accurate deep potential. The use of 32 B atoms to generate the training set is a choice that balances efficiency and accuracy.

In the DPMD simulations, we adopt the NVT ensemble with the Nosé–Hoover thermostat.69,70 We use the LAMMPS package.71 The number of B atoms ranges from 32 to 16 384. Time steps of 0.07 and 0.01 fs are set for the systems at 86 and 350 eV, respectively. We perform 400 000 steps of DPMD simulations to yield the radial distribution functions, the static structure factors, and the dynamic structure factors. Furthermore, 10⁶ molecular dynamics steps are performed to compute the shear viscosity.

After the BOMD trajectories are generated, the radial distribution function can be evaluated according to $g (r) = \frac{V}{4 π r^{2} N (N - 1)} ⟨\sum_{i = 1}^{N} \sum_{j = 1, j \neq i}^{N} δ (r - | r_{i} - r_{j} |)⟩,$ (16)where V is the cell volume, N is the number of atoms, r_i and r_j are the atomic coordinates of atoms i and j, and ⟨⋯⟩ denotes the ensemble average.

We plot g(r) of warm dense B with a density of 2.46 g/cm³ 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 the BOMD simulations.

Radial distribution functions g(r) of B systems with a density of 2.46 g/cm3 at temperatures of (a) 86 eV and (b) 350 eV. The g(r) obtained by Ext-FPMD with the local density approximation (LDA) functional comes from Blanchet et al.19 The SDFT calculations are performed with the PBE59 and corrKSDT65 XC functionals. The number of B atoms is set to 32 in the first-principles molecular dynamics simulations. DPMD denotes the model trained by the traditional DP method50 and DPMD-T the model trained by the TDDP method.51N is the number of B atoms in a cell, which ranges from 32 to 16 384 in the deep-potential-based simulations.

Fig. 6. Radial distribution functions g(r) of B systems with a density of 2.46 g/cm³ at temperatures of (a) 86 eV and (b) 350 eV. The g(r) obtained by Ext-FPMD with the local density approximation (LDA) functional comes from Blanchet et al.¹⁹ The SDFT calculations are performed with the PBE⁵⁹ and corrKSDT⁶⁵ XC functionals. The number of B atoms is set to 32 in the first-principles molecular dynamics simulations. DPMD denotes the model trained by the traditional DP method⁵⁰ and DPMD-T the model trained by the TDDP method.⁵¹N is the number of B atoms in a cell, which ranges from 32 to 16 384 in the deep-potential-based simulations.

下载图片查看所有图片

The results reveal the following.

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 XC functional59 and the finite-temperature XC functional corrKSDT,65 which indicates that temperature effects in the XC functional are not significant. Third, as expected, the g(r) is not smooth, owing to the limited number of molecular dynamics steps (4000). However, by employing the DPMD method, we obtain not only an accurate g(r) that agrees well with the SDFT results but also a smooth g(r), since a larger number of atoms (108 to 16 384) and a longer trajectory (400 000 steps) are considered. Importantly, size effects can be largely mitigated, as evidenced by the convergence of g(r) at around 108 atoms.

3.5 E. Ion–ion static structure factors

The ion–ion static structure factor S(q) measured from diffraction experiments72,73 contains information regarding the spatial arrangement of particles in a material. The formula for S(q) is $S (q) = \frac{1}{N} ⟨\sum_{i = 1}^{N} \sum_{j = 1}^{N} e^{i q \cdot (r_{i} - r_{j})}⟩,$ (17)where N is the number of atoms, i and j denote atoms, and q is the wave vector.

Here, we perform BOMD simulations on a 32-atom cell by SDFT with the PBE59 and corrKSDT65 XC functionals. Moreover, we employ DPMD to calculate S(q) for cells containing 32, 108, 256, 2048, and 16 384 B atoms with a density of 2.46 g/cm³. The results for systems at 86 and 350 eV are illustrated in Figs. 7(a) and 7(b), respectively.

Static structure factors S(q) of B with a density of 2.46 g/cm3 at temperatures of (a) 86 eV and (b) 350 eV. The SDFT calculations are performed with the PBE59 and corrKSDT65 XC functionals. The number of B atoms is set to 32 in the first-principles molecular dynamics simulations. DPMD denotes the model trained by the traditional DP method50 and DPMD-T the model trained by the TDDP method.51N is the number of B atoms in a cell and ranges from 32 to 16 384 in the deep-potential-based simulations.

Fig. 7. Static structure factors S(q) of B with a density of 2.46 g/cm³ at temperatures of (a) 86 eV and (b) 350 eV. The SDFT calculations are performed with the PBE⁵⁹ and corrKSDT⁶⁵ XC functionals. The number of B atoms is set to 32 in the first-principles molecular dynamics simulations. DPMD denotes the model trained by the traditional DP method⁵⁰ and DPMD-T the model trained by the TDDP method.⁵¹N is the number of B atoms in a cell and ranges from 32 to 16 384 in the deep-potential-based simulations.

下载图片查看所有图片

It is noteworthy that the data points of S(q) generated by SDFT exhibit oscillations due to the limited number of molecular dynamics steps (4000). However, the DPMD simulations offer better converged results, since 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 about S(q), which signifies the long-ranged order of a system.

3.6 F. Ion–ion dynamic structure factors

The collective dynamics of ion density fluctuations can be characterized by the ion–ion dynamic structure factor S(q, ω), which is experimentally measurable74 and plays a crucial role in investigating ion dynamics, including collective modes,75 dissipation processes,76 and others. In practice, S(q, ω) can be computed from the intermediate scattering function F(q, t) via Fourier transformation using the following formula: $S (q, ω) = \frac{1}{2 π} \int_{- \infty}^{+ \infty} F (q, t) e^{i ω t} d t .$ (18)Here F(q, t) takes the form $F (q, t) = \frac{1}{N} ⟨ ρ (q, t) ρ (- q, 0) ⟩,$ (19)where N is the number of ions, and ρ(q, t) is defined as $ρ (q, t) = \sum_{i = 1}^{N} e^{i q \cdot r_{i} (t)},$ (20)where r_i(t) denotes the atomic coordinates for atom i at time t.

Figures 8(a) and 8(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) and three sizes of cells are tested (256, 2048, and 16 384 atoms). We find that 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.

Intermediate scattering functions F(q, t) of B with a density of 2.46 g/cm3 as calculated from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted in DPMD simulations at two temperatures of (a) 86 eV and (b) 350 eV. Three wave vectors are chosen: q = 0.51, 1.02, and 2.50 Å−1.

Fig. 8. Intermediate scattering functions F(q, t) of B with a density of 2.46 g/cm³ as calculated from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted in DPMD simulations at two temperatures of (a) 86 eV and (b) 350 eV. Three wave vectors are chosen: q = 0.51, 1.02, and 2.50 Å⁻¹.

下载图片查看所有图片

Next, we obtain the ion–ion dynamic structure factors S(q, ω) of warm dense B by performing the Fourier transform of F(q, t). The results associated with wave vectors q = 0.51, 1.02, and 2.50 Å⁻¹ are shown in Figs. 9(a)–9(c), respectively. In each figure, two temperatures (86 and 350 eV) and three system sizes (256, 2048, and 16 384 atoms) are adopted.

Ion–ion dynamic structure factors S(q, ω) of B with a density of 2.46 g/cm3 as computed from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted. The wave vectors q are chosen to be (a) 0.51 Å−1, (b) 1.02 Å−1, and (c) 2.50 Å−1.

Fig. 9. Ion–ion dynamic structure factors S(q, ω) of B with a density of 2.46 g/cm³ as computed from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted. The wave vectors q are chosen to be (a) 0.51 Å⁻¹, (b) 1.02 Å⁻¹, and (c) 2.50 Å⁻¹.

下载图片查看所有图片

For the wave vector, q = 0.51 Å⁻¹, Fig. 9(a) shows well-pronounced ion-acoustic modes located at ω = 206.78 for 86 and 486.80 meV for 350 eV. When q increases to 1.02 Å⁻¹ in Fig. 9(b), the peak of S(q, ω) becomes lower, and its location shifts to 324.95 for 86 and 616.53 meV for 350 eV. Notably, the ion-acoustic mode S(q, ω) disappears when q = 2.50 Å⁻¹, because the noncollective mode dominates at large q. These results for S(q, ω) demonstrate that the DP method can predict the long-ranged structural and time correlation of WDM. For high temperatures up to hundreds of eV, there are experimental measurements of ion–ion static and dynamic structure factors via x-ray Thomson scattering for materials such as CH77 and Be,78 but, to the best of our knowledge, no experimental data are available for B at temperatures of 86 and 350 eV.

3.7 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,79,80 according to which the shear viscosity is linked to the integral of the stress auto-correlation function in the form $η = \frac{V}{3 k_{B} T} \int_{0}^{+ \infty} [\sum_{α β} ⟨ P_{α β} (0) P_{α β} (t) ⟩] d t,$ (21)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 autocorrelation functions of B at a density of 2.46 g/cm³ and temperatures of 86 and 350 eV are displayed in Figs. 10(a) and 10(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 16 384.

Stress autocorrelation functions [Eq. (21)] of warm dense B with a density of 2.46 g/cm3 at temperatures of (a) 86 eV and (b) 350 eV. (c) Shear viscosity of B. DPMD simulations are used, with cells containing 32, 108, 256, 864, 2048, 6912, and 16 384 atoms. Error bars represent standard deviations.

Fig. 10. Stress autocorrelation functions [Eq. (21)] of warm dense B with a density of 2.46 g/cm³ at temperatures of (a) 86 eV and (b) 350 eV. (c) Shear viscosity of B. DPMD simulations are used, with cells containing 32, 108, 256, 864, 2048, 6912, and 16 384 atoms. Error bars represent standard deviations.

下载图片查看所有图片

During DPMD simulations, each system is first relaxed for 50 000 steps. Next, 10⁶ steps of molecular dynamics simulations are performed to calculate the stress auto-correlation function. In detail, the trajectory length is 70 for 86 and 10 ps for 350 eV. We take values from 0.105 to 0.305 and from 0.05 to 0.1 ps to compute the averaged shear viscosity for the systems at 86 and 350 eV, respectively, 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 exhibit no substantial finite-size effects for the shear viscosity of B, which is consistent with previous conclusions.81–83

There are other formulas that can be used to predict the shear viscosity of materials.

For example, we note that a recent work proposes an extended random-walk shielding-potential viscosity model (ext-RWSP-VM)84,85 to evaluate the shear viscosity of materials in WDM and HDP states. The viscosity is evaluated as $η = \frac{\sqrt{3 m k_{B} T}}{π d^{4}} I,$ (22)where d is the collision diameter according to the hard-sphere 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 at temperatures of 86 and 350 eV, respectively. Further details can be found in the supplementary material. Interestingly, the data are close to our first-principles results.

In addition, we find that the shear viscosity of a plasma can also be described by the approximate formula86 $η = \frac{2}{3 \sqrt{π}} \frac{\sqrt{m k_{B} T}}{σ_{o}},$ (23)where m is the atomic mass and σ_o is the total collision cross section (∼10⁻²⁰ m²). Thus, the estimated shear viscosities of B are ∼18.7 and 37.8 mPa · s at 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.86 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 shear viscosities calculated from DPMD are also consistent with the approximate values obtained using Eq. (23).

4 IV. CONCLUSIONS

Simulating WDM with first-principles accuracy has long been challenging, owing to the partial occupation of a large number of high-energy KS eigenstates and the resulting limitations on the time and space scales. Our work has suggested that the advent of the SDFT method and machine-learning-based molecular dynamics could be of great help in overcoming the difficulties. The SDFT method described in this work has been implemented with a plane-wave basis set and the k-point sampling method, which has been enabled in the ABACUS package. In this work, we have 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 have substantially extended the time and space scales for simulating warm dense B and have also reduced the finite size effect. In addition, we have 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 has validated the combination of stochastic density functional theory with machine learning techniques to study high-temperature systems. It has also provided new insights into the properties of WDM. 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.

5 SUPPLEMENTARY MATERIAL

6 ACKNOWLEDGMENTS

Acknowledgment. We thank Xinyu Zhang for helpful discussions. We thank Hang Zhang for proofreading the manuscript. We thank the electronic structure team (from AI for Science Institute, Beijing) for improving the ABACUS package from various aspects. This work is supported by the National Natural Science Foundation of China under Grant Nos. 12122401 and 12074007. The numerical simulations were performed on the High-Performance Computing Platform of the Beijing Super Cloud Computing Center and the Bohrium cloud computing platform of DP Technology Co., Ltd.

References

[1] T.Guillot. Interiors of giant planets inside and outside the solar system. Science, 1999, 286: 72-77.

[2] J.Bates, R.Betti, E.Campbell, T.Chapman, T.Collins, D.Froula, V.Goncharov, C.Goyon, M.Hohenberger, I.Igumenshchev, M.Karasik, J.Knauer, J.Marozas, F.Marshall, L.Masse, A.Maximov, R.McCrory, P.Michel, J.Myatt, K.Obenschain, S.Obenschain, J.Oh, P.Radha, S.Regan, S.Reyes, M.Rosenberg, S.Ross, T.Sangster, A.Schmitt, A.Sefkow, W.Seka, A.Shvydky, A.Solodov, D.Turnbull, B.Van Wonterghem, J.Weaver. Laser-direct-drive program: Promise, challenge, and path forward. Matter Radiat. Extremes, 2017, 2: 37-54.

[3] E. I.Moses. Ignition on the National Ignition Facility: A path towards inertial fusion energy. Nucl. Fusion, 2009, 49: 104022.

[4] F.Graziani, M. P.Desjarlais, R.Redmer, and S. B.Trickey, Frontiers and Challenges in Warm Dense Matter (SpringerCham, 2014).

[5] R. P.Drake, B. A.Remington, D. D.Ryutov. Experimental astrophysics with high power lasers and Z pinches. Rev. Mod. Phys., 2006, 78: 755-807.

[6] P.Hohenberg, W.Kohn. Inhomogeneous electron gas. Phys. Rev., 1964, 136: B864-B871.

[7] W.Kohn, L. J.Sham. Self-consistent equations including exchange and correlation effects. Phys. Rev., 1965, 140: A1133-A1138.

[8] D. M.Ceperley, B.Militzer. Path integral Monte Carlo calculation of the deuterium hugoniot. Phys. Rev. Lett., 2000, 85: 1890-1893.

[9] D. M.Ceperley, L. A.Collins, J. D.Johnson, J. D.Kress, S.Mazevet, B.Militzer. Calculation of a deuterium double shock Hugoniot from ab initio simulations. Phys. Rev. Lett., 2001, 87: 275502.

[10] V. N.Goncharov, S. X.Hu, B.Militzer, S.Skupsky. Strong coupling and degeneracy effects in inertial confinement fusion implosions. Phys. Rev. Lett., 2010, 104: 235003.

[11] K. P.Driver, B.Militzer. All-electron path integral Monte Carlo simulations of warm dense matter: Application to water and carbon plasmas. Phys. Rev. Lett., 2012, 108: 115502.

[12] X.-T.He, C.Wang, P.Zhang. Ab initio simulations of dense helium plasmas. Phys. Rev. Lett., 2011, 106: 145002.

[13] B. J. B.Crowley, G.Gregori, J. W. O.Harris, L. K.Pattison, S.Richardson, T. G.White. Orbital-free density-functional theory simulations of the dynamic structure factor of warm dense aluminum. Phys. Rev. Lett., 2013, 111: 175002.

[14] D.Kang, K.Luo, K.Runge, S. B.Trickey. Two-temperature warm dense hydrogen as a test of quantum protons driven by orbital-free density functional theory electronic forces. Matter Radiat. Extremes, 2020, 5: 064403.

[15] M.Chen, Q.Liu, D.Lu. Structure and dynamics of warm dense aluminum: A molecular dynamics study with density functional theory and deep potential. J. Phys.: Condens. Matter, 2020, 32: 144002.

[16] X. T.He, W.Kang, H.Wang, P.Zhang, S.Zhang. Extended application of Kohn-Sham first-principles molecular dynamics method with plane wave approximation at high energy—From cold materials to hot dense plasmas. Phys. Plasmas, 2016, 23: 042707.

[17] C.Gao, X. T.He, W.Kang, D.Li, X.Liu, C.Wang, P.Zhang, S.Zhang, W.Zhang, X.Zhang. Equations of state of poly-α-methylstyrene and polystyrene: First-principles calculations versus precision measurements. Phys. Rev. B, 2021, 103: 174111.

[18] A.Blanchet, J.Clérouin, F.Soubiran, M.Torrent. Extended first-principles molecular dynamics model for high temperature simulations in the ABINIT code: Application to warm dense aluminum. Comput. Phys. Commun., 2022, 271: 108215.

[19] A.Blanchet, J.Clérouin, F.Soubiran, M.Torrent. Extended first-principles molecular dynamics simulations of hot dense boron: Equation of state and ionization. Contrib. Plasma Phys., 2022, 62: e202100234.

[20] R.Baer, D.Neuhauser, E.Rabani. Self-averaging stochastic Kohn-Sham density-functional theory. Phys. Rev. Lett., 2013, 111: 106402.

[21] R.Baer, Y.Cytter, D.Neuhauser, E.Rabani. Stochastic density functional theory at finite temperatures. Phys. Rev. B, 2018, 97: 115207.

[22] R.Baer, M. D.Fabian, D.Neuhauser, E.Rabani, B.Shpiro. Stochastic density functional theory. Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9: e1412.

[23] R.Baer, D.Neuhauser, E.Rabani. Stochastic vector techniques in ground-state electronic structure. Annu. Rev. Phys. Chem., 2022, 73: 255-272.

[24] L. A.Collins, V.Sharma, A. J.White. Stochastic and mixed density functional theory within the projector augmented wave formalism for simulation of warm dense matter. Phys. Rev. E, 2023, 108: L023201.

[25] N. D.Mermin. Thermal properties of the inhomogeneous electron gas. Phys. Rev., 1965, 137: A1441-A1443.

[26] T. W.Barbee, M. P.Surh, L. H.Yang. First principles molecular dynamics of dense plasmas. Phys. Rev. Lett., 2001, 86: 5958-5961.

[27] C.Wang, P.Zhang. Wide range equation of state for fluid hydrogen from density functional theory. Phys. Plasmas, 2013, 20: 092703.

[28] J.Daligault, T.Sjostrom. Fast and accurate quantum molecular dynamics of dense plasmas across temperature regimes. Phys. Rev. Lett., 2014, 113: 155006.

[29] A.Cangi, T.Dornheim, L.Fiedler, K.Jiang, Z. A.Moldabekov, M.Pavanello, X.Shao. Accelerating equilibration in first-principles molecular dynamics with orbital-free density functional theory. Phys. Rev. Res., 2022, 4: 043033.

[30] K. P.Driver, F.González-Cataldo, B.Militzer, F. m. c.Soubiran, S.Zhang. First-principles equation of state database for warm dense matter computation. Phys. Rev. E, 2021, 103: 013203.

[31] T.Dornheim. Fermion sign problem in path integral Monte Carlo simulations: Quantum dots, ultracold atoms, and warm dense matter. Phys. Rev. E, 2019, 100: 023307.

[32] V. V.Karasiev, K.Luo, S. B.Trickey. A simple generalized gradient approximation for the noninteracting kinetic energy density functional. Phys. Rev. B, 2018, 98: 041111.

[33] L. A.Collins, A. J.White. Fast and universal Kohn-Sham density functional theory algorithm for warm dense matter to hot dense plasma. Phys. Rev. Lett., 2020, 125: 055002.

[34] M.Chen, G.-C.Guo, L.He. Systematically improvable optimized atomic basis sets for ab initio calculations. J. Phys.: Condens. Matter, 2010, 22: 445501.

[35] M.Chen, L.He, P.Li, L.Lin, P.Lin, X.Liu, X.Ren, C.Yang. Large-scale ab initio simulations based on systematically improvable atomic basis. Comput. Mater. Sci., 2016, 112: 503-517.

[36] M.Chen, Q.Liu. Plane-wave-based stochastic-deterministic density functional theory for extended systems. Phys. Rev. B, 2022, 106: 125132.

[37] J.Behler. Neural network potential-energy surfaces in chemistry: A tool for large-scale simulations. Phys. Chem. Chem. Phys., 2011, 13: 17930-17955.

[38] J.Behler, T.Morawietz. A density-functional theory-based neural network potential for water clusters including van der waals corrections. J. Phys. Chem. A, 2013, 117: 7356-7366.

[39] A. P.Bartók, G.Csányi. Gaussian approximation potentials: A brief tutorial introduction. Int. J. Quantum Chem., 2015, 115: 1051-1057.

[40] R.Car, R. A.DiStasio Jr, H.-Y.Ko, B.Santra, H.Wang, L.Zhang. Isotope effects in liquid water via deep potential molecular dynamics. Mol. Phys., 2019, 117: 3269-3281.

[41] R.Batra, Y.Mishin, G. P. P.Pun, R.Ramprasad. Physically informed artificial neural networks for atomistic modeling of materials. Nat. Commun., 2019, 10: 2339.

[42] K.Barros, C.Devereux, O.Isayev, N.Lubbers, B. T.Nebgen, A. E.Roitberg, J. S.Smith, S.Tretiak, R.Zubatyuk. Approaching coupled cluster accuracy with a general-purpose neural network potential through transfer learning. Nat. Commun., 2019, 10: 2903.

[43] B.Chen, J.Dai, D.Kang, H.Wang, X.Yu, Q.Zeng, S.Zhang. Towards large-scale and spatiotemporally resolved diagnosis of electronic density of states by deep learning. Phys. Rev. B, 2022, 105: 174109.

[44] R.Car, P. G.Debenedetti, T. E.Gartner, A. Z.Panagiotopoulos, P. M.Piaggi, L.Zhang. Signatures of a liquid–liquid transition in an ab initio deep neural network model for water. Proc. Natl. Acad. Sci. U. S. A., 2020, 117: 26040-26046.

[45] M.French, R.Redmer, H. R.Rüter, M.Sch?rner. Extending ab initio simulations for the ion-ion structure factor of warm dense aluminum to the hydrodynamic limit using neural network potentials. Phys. Rev. B, 2022, 105: 174310.

[46] A. D.Baczewski, A.Cangi, R.Redmer, M.Sch?rner, B. B. L.Witte. Ab initio study of shock-compressed copper. Phys. Rev. B, 2022, 106: 054304.

[47] M.Bussmann, A.Cangi, L.Fiedler, K.Shah. Deep dive into machine learning density functional theory for materials science and chemistry. Phys. Rev. Mater., 2022, 6: 040301.

[48] A.Cangi, S.Kumar, M.Lokamani, S.Nikolov, K.Ramakrishna, H.Tahmasbi, J.Tranchida, M. A.Wood. Transferable interatomic potential for aluminum from ambient conditions to warm dense matter. Phys. Rev. Res., 2023, 5: 033162.

[49] R.Car, W.E, J.Han, H.Wang, L.Zhang. Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics. Phys. Rev. Lett., 2018, 120: 143001.

[50] W.E, J.Han, H.Wang, L.Zhang. DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics. Comput. Phys. Commun., 2018, 228: 178-184.

[51] M.Chen, C.Gao, Q.Liu, H.Wang, L.Zhang, Y.Zhang. Warm dense matter simulation via electron temperature dependent deep potential molecular dynamics. Phys. Plasmas, 2020, 27: 122704.

[52] W.E, D. J.Srolovitz, H.Wang, T.Wen, L.Zhang. Deep potentials for materials science. Mater. Futures, 2022, 1: 022601.

[53] R.Car, M.Chen, W.E, W.Jia, L.Lin, D.Lu, H.Wang, L.Zhang. 86 PFLOPS deep potential molecular dynamics simulation of 100 million atoms with ab initio accuracy. Comput. Phys. Commun., 2021, 259: 107624.

[54] W.Jia, H.Wang, M.Chen, D.Lu, L.Lin, R.Car, W.E, and L.Zhang, Pushing the Limit of Molecular Dynamics with Ab Initio Accuracy to 100 Million Atoms with Machine Learning (IEEE Press, 2020).

[55] B.Chen, J.Dai, T.Gao, D.Kang, H.Wang, Y.Yao, X.Yu, Q.Zeng, S.Zhang. Ab initio validation on the connection between atomistic and hydrodynamic description to unravel the ion dynamics of warm dense matter. Phys. Rev. Res., 2021, 3: 033116.

[56] M.Chen, J.Li, Q.Liu. Thermal transport by electrons and ions in warm dense aluminum: A combined density functional theory and deep potential study. Matter Radiat. Extremes, 2021, 6: 026902.

[57] R.Baer, M.Head-Gordon. Sparsity of the density matrix in Kohn-Sham density functional theory and an assessment of linear system-size scaling methods. Phys. Rev. Lett., 1997, 79: 3962-3965.

[58] P.Kingma and J.Ba, “Adam: A method for stochastic optimization,” [cs.LG] (2017).

[59] K.Burke, M.Ernzerhof, J. P.Perdew. Generalized gradient approximation made simple. Phys. Rev. Lett., 1996, 77: 3865-3868.

[60] M.Chen, Y.Chen, W.Jia, W.Jiang, D.Lu, H.Wang, L.Zhang. DP compress: A model compression scheme for generating efficient deep potential models. J. Chem. Theory Comput., 2022, 18: 5559-5567.

[61] D. R.Hamann. Optimized norm-conserving Vanderbilt pseudopotentials. Phys. Rev. B, 2013, 88: 085117.

[62] F.Gygi, M.Schlipf. Optimization algorithm for the generation of ONCV pseudopotentials. Comput. Phys. Commun., 2015, 196: 36-44.

[63] C.Gao, X. T.He, W.Kang, C.Wang, P.Zhang, S.Zhang. Validity boundary of orbital-free molecular dynamics method corresponding to thermal ionization of shell structure. Phys. Rev. B, 2016, 94: 205115.

[64] K.Caspersen, P. M.Celliers, D.Erskine, J.Gaffney, M. C.Gregor, A.Lazicki, R. A.London, B.Militzer, J.Nilsen, T.Ogitsu, P. A.Sterne, D.Swift, H. D.Whitley, L. H.Yang, S.Zhang. Theoretical and experimental investigation of the equation of state of boron plasmas. Phys. Rev. E, 2018, 98: 023205.

[65] J. W.Dufty, V. V.Karasiev, S. B.Trickey. Nonempirical semilocal free-energy density functional for matter under extreme conditions. Phys. Rev. Lett., 2018, 120: 076401.

[66] R.Car, W.E, D.-Y.Lin, H.Wang, L.Zhang. Active learning of uniformly accurate interatomic potentials for materials simulation. Phys. Rev. Mater., 2019, 3: 023804.

[67] L.Cao, M.Xu, J.Zeng, J. Z. H.Zhang, T.Zhu. Complex reaction processes in combustion unraveled by neural network-based molecular dynamics simulation. Nat. Commun., 2020, 11: 5713.

[68] M.Chen, T.Chen, H.Geng, J.Liu, H.Wang, F.Yuan, L.Zhang. Modeling the high-pressure solid and liquid phases of tin from deep potentials with ab initio accuracy. Phys. Rev. Mater., 2023, 7: 053603.

[69] S.Nosé. A unified formulation of the constant temperature molecular dynamics methods. J. Chem. Phys., 1984, 81: 511-519.

[70] W. G.Hoover. Canonical dynamics: Equilibrium phase-space distributions. Phys. Rev. A, 1985, 31: 1695-1697.

[71] S.Plimpton. Fast parallel algorithms for short-range molecular dynamics. J. Comput. Phys., 1995, 117: 1-19.

[72] A. J.Greenfield, J.Wellendorf, N.Wiser. X-ray determination of the static structure factor of liquid Na and K. Phys. Rev. A, 1971, 4: 1607-1616.

[73] P.Martel, V. F.Sears, E. C.Svensson, A. D. B.Woods. Neutron-diffraction study of the static structure factor and pair correlations in liquid ⁴He. Phys. Rev. B, 1980, 21: 3638-3651.

[74] H.K?hlert. Dynamic structure factor of strongly coupled Yukawa plasmas with dissipation. Phys. Plasmas, 2019, 26: 063703.

[75] V. M.Giordano, G.Monaco. Fingerprints of order and disorder on the high-frequency dynamics of liquids. Proc. Natl. Acad. Sci. U. S. A., 2010, 107: 21985-21989.

[76] L. B.Fletcher, D. O.Gericke, S. H.Glenzer, G.Gregori, N. J.Hartley, P.Mabey, S.Richardson, J.Vorberger, T. G.White. A strong diffusive ion mode in dense ionized matter predicted by Langevin dynamics. Nat. Commun., 2017, 8: 14125.

[77] B.Bachmann, R. A.Baggott, D. A.Chapman, G. W.Collins, T.D?ppner, R. W.Falcone, D. O.Gericke, S. H.Glenzer, J. A.Hawreliak, D. H.Kalantar, D.Kraus, A. L.Kritcher, O. L.Landen, S.Le Pape, T.Ma, P.Neumayer, J.Nilsen, D. C.Swift. X-ray scattering measurements on imploding CH spheres at the National Ignition Facility. Phys. Rev. E, 2016, 94: 011202.

[78] B.Bachmann, R. A.Baggott, M.Bethkenhagen, M. P.B?hme, D. A.Chapman, L.Divol, T.D?ppner, R. W.Falcone, L. B.Fletcher, D. O.Gericke, S. H.Glenzer, D.Kraus, O. L.Landen, M. J.MacDonald, P.Neumayer, R.Redmer, A. M.Saunders, M.Sch?rner, P. A.Sterne, J.Vorberger, B. B. L.Witte, A.Yi. Observing the onset of pressure-driven K-shell delocalization. Nature, 2023, 618: 270-275.

[79] M. S.Green. Markoff random processes and the statistical mechanics of time-dependent phenomena. II. Irreversible processes in fluids. J. Chem. Phys., 1954, 22: 398-413.

[80] R.Kubo. Statistical-mechanical theory of irreversible processes. I. General theory and simple applications to magnetic and conduction problems. J. Phys. Soc. Jpn., 1957, 12: 570-586.

[81] E. A.Carter, M.Chen, P. G.Debenedetti, A. Z.Panagiotopoulos, F. H.Stillinger, J. R.Vella. Liquid li structure and dynamics: A comparison between OFDFT and second nearest-neighbor embedded-atom method. AIChE J., 61: 2841-2853.

[82] E. A.Carter, M.Chen, P. G.Debenedetti, A. Z.Panagiotopoulos, F. H.Stillinger, J. R.Vella. Structural and dynamic properties of liquid tin from a new modified embedded-atom method force field. Phys. Rev. B, 2017, 95: 064202.

[83] S.Baroni, R.Car, C.Malosso, D.Tisi, L.Zhang. Viscosity in water from first-principles and deep-neural-network simulations. npj Comput. Mater., 2022, 8: 139.

[84] Y.Cheng, X.Gao, Y.Hou, Q.Li, H.Liu, Y.Liu, X.Meng, H.Song, J.Wang, J.Yuan. Random-walk shielding-potential viscosity model for warm dense metals. Phys. Rev. E, 2022, 106: 014142.

[85] Y.Cheng, X.Gao, Q.Li, Y.Liu, H.Song, and H.Liu, “Extended application of random-walk shielding-potential viscosity model of metals in wide temperature region,” [cond-mat.stat-mech] (2023).

[86] E. P. M. I.Boulos and P. L.Fauchais, Handbook of Thermal Plasmas (Springer, Cham, 2023).

1 I. INTRODUCTION

2 II. COMPUTATIONAL METHODS

2.1 A. Stochastic and mixed stochastic–deterministic density functional theory

2.2 B. Deep potential molecular dynamics

3 III. RESULTS AND DISCUSSION

3.1 A. Statistical errors of SDFT

3.2 B. Comparison of SDFT and KSDFT results

3.3 C. High-temperature calculations by SDFT

3.4 D. Radial distribution functions

3.5 E. Ion–ion static structure factors

3.6 F. Ion–ion dynamic structure factors

3.7 G. Shear viscosities

4 IV. CONCLUSIONS

5 SUPPLEMENTARY MATERIAL

6 ACKNOWLEDGMENTS

Tao Chen, Qianrui Liu, Yu Liu, Liang Sun, Mohan Chen. Combining stochastic density functional theory with deep potential molecular dynamics to study warm dense matter[J]. Matter and Radiation at Extremes, 2024, 9(1): 015604.

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

1 I. INTRODUCTION

2 II. COMPUTATIONAL METHODS

2.1 A. Stochastic and mixed stochastic–deterministic density functional theory

2.2 B. Deep potential molecular dynamics

3 III. RESULTS AND DISCUSSION

3.1 A. Statistical errors of SDFT

3.2 B. Comparison of SDFT and KSDFT results

3.3 C. High-temperature calculations by SDFT

3.4 D. Radial distribution functions

3.5 E. Ion–ion static structure factors

3.6 F. Ion–ion dynamic structure factors

Fig. 9. Ion–ion dynamic structure factors S(q, ω) of B with a density of 2.46 g/cm³ as computed from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted. The wave vectors q are chosen to be (a) 0.51 Å⁻¹, (b) 1.02 Å⁻¹, and (c) 2.50 Å⁻¹.

3.7 G. Shear viscosities

4 IV. CONCLUSIONS

5 SUPPLEMENTARY MATERIAL

6 ACKNOWLEDGMENTS

Article Outline

关于本站 Cookie 的使用提示

全站搜索

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

1 I. INTRODUCTION

2 II. COMPUTATIONAL METHODS

2.1 A. Stochastic and mixed stochastic–deterministic density functional theory

2.2 B. Deep potential molecular dynamics

3 III. RESULTS AND DISCUSSION

3.1 A. Statistical errors of SDFT

3.2 B. Comparison of SDFT and KSDFT results

3.3 C. High-temperature calculations by SDFT

3.4 D. Radial distribution functions

3.5 E. Ion–ion static structure factors

3.6 F. Ion–ion dynamic structure factors

Fig. 9. Ion–ion dynamic structure factors S(q, ω) of B with a density of 2.46 g/cm3 as computed from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted. The wave vectors q are chosen to be (a) 0.51 Å−1, (b) 1.02 Å−1, and (c) 2.50 Å−1.

3.7 G. Shear viscosities

4 IV. CONCLUSIONS

5 SUPPLEMENTARY MATERIAL

6 ACKNOWLEDGMENTS

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

Fig. 9. Ion–ion dynamic structure factors S(q, ω) of B with a density of 2.46 g/cm³ as computed from DPMD trajectories. Three system sizes (256, 2048, and 16 384 atoms) are adopted. The wave vectors q are chosen to be (a) 0.51 Å⁻¹, (b) 1.02 Å⁻¹, and (c) 2.50 Å⁻¹.