Transport coefficients of warm dense matter from Kohn-Sham density functional theory Open Access

The warm dense matter (WDM) regime exists at the intersection of plasma and condensed matter physics, with solid densities $ρ = 10 − 2 − 10 4$ g/cm³ and relatively high temperatures $T = 1 − 100$ eV.¹ Accurate material models are of fundamental importance for many areas of research in planetary science,^2,3 astrophysics,⁴ inertial confinement fusion,^5–7 among others. Yet, despite significant efforts at experimental high energy density (HED) facilities such as the National Ignition Facility (NIF) at Lawrence Livermore National Laboratory (LLNL)⁸ or the Z-machine at Sandia National Laboratories (SNL),⁹ accurate measurements of transport coefficients in the WDM regime remain exceptionally difficult and therefore are relatively rare. As a result, there is a strong reliance on theoretical approaches for equations of state and transport properties.

From a theoretical point of view, the WDM regime is challenging due to the importance of both electron degeneracy $Θ = k B T e E F < 1$ and strong ion coupling $Γ = ( Z * ) 2 e 2 / a ws T i ≫ 1$⁠, where k_B is the Boltzmann constant, T_e is the electronic temperature, E_F is the Fermi energy, $Z *$ is the ionization of the ions, $a ws$ is the Wigner–Seitz radius, and T_i is the ionic temperature.¹ Models of plasma that are accurate at extremely high temperature often neglect electron correlation, atomic shell structure, and other quantum phenomena that are essential in the WDM regime.^10,11 In contrast, commonly used ab initio methods like Kohn–Sham density functional theory based molecular dynamics (DFT-MD) explicitly treat many-body ion dynamics and electronic structure. Typical approaches rely on approximate semi-local exchange–correlation functionals—such as the localized density approximations (LDA) using the Perdew–Zunger¹² or the Vosko–Wilk–Nussair¹³ parameterization of the Ceperley–Alder¹⁴ homogeneous electron gas data or generalized gradient approximations such as Perdew–Burke–Ernzerhof (PBE)^15,16 models—that capture some electron correlation effects. The reliance on the Kohn–Sham single particle description¹⁷ requires a large number of orbitals with small occupation at high temperatures due to the Fermi–Dirac distribution and therefore scales poorly with temperature.¹⁸ We note that this scaling can be partially mitigated by utilizing a stochastic orbital representation or mixed basis sets instead of a pure Kohn–Sham scheme.^19,20 Additionally, the use of ground state exchange–correlation functionals may not be adequate for materials in the WDM regime and can be improved through the use of finite-temperature functionals;^21–23 however, these are not yet widely available in ab initio codes such as Vienna ab initio Simulation Package (VASP),^24–26 which we use here. Nevertheless, where the calculations are feasible, DFT-MD is considered a benchmark method in the HED regime for systems where the ions and electrons are in equilibrium.

Obtaining accurate transport properties from a DFT-MD simulation requires many technical considerations. Well-converged estimates for a particular property require understanding its sensitivity to many details of the simulation including, e.g., time step, energy cutoff, system size, etc. Moreover, sensitivities to these choices can be vastly different based on the property of interest. For example, the total energy and pressure exhibit relatively fast convergence with respect to the number of bands, k-points, the length of simulation, the number of atoms, and even the choice of thermostat. In contrast, the DC limit of the electrical conductivity exhibits much greater sensitivity to most of these variables.^27,28 Accurate estimates of the shear viscosity require long simulations with the appropriate thermostat and appreciation of its unfavorable statistical properties.^29,30 While the above considerations are generally well known in principle, in practice it can be difficult to reproduce results from the literature because such details are often briefly mentioned or omitted entirely.

In this paper, we present our contributions for the “priority 1” cases for the Second Charged-Particle Transport Coefficient Code Comparison Workshop (TCCW2) using ab initio Molecular dynamics (MD) based on Kohn–Sham Density Functional Theory (DFT), summarized in Table I. In the first workshop, various methods were applied to H, C, and CH for a variety of densities and temperatures in the WDM regime. For the second workshop, a priority system was introduced to help guide contributors about what systems were most important. The priority 1 systems were chosen to enable comparisons between the first and second workshops, as well as introduce more complicated systems at densities and temperatures such that DFT-MD was still applicable; this allows DFT-MD methods to be used to help validate plasma and kinetic models, which can be extended to much more extreme conditions. We describe in detail how to calculate the various transport coefficients including the electronic contributions to the DC electrical conductivity and thermal conductivity, as well as the shear viscosity within the DFT-MD framework. In addition, we discuss optical properties such as the index of refraction, absorption coefficient, and reflectivity, which can be obtained from the electrical conductivity. In Sec. II, we discuss the theoretical considerations for each quantity as well as the detailed approach used to estimate these properties. In Sec. III, results are presented for each of the systems outlined in the priority 1 cases of the TCCW2.

Producing robust and accurate estimates of transport properties as well as their corresponding statistical uncertainties from DFT-MD data requires consideration of several aspects of the calculation and underlying theory. In this section, we discuss in detail both the theory and practical convergence for a variety of transport and optical properties. Our goal is to clearly demonstrate how we have chosen to analyze and converge each quantity in order to improve reproducibility.

In Sec. II A, we discuss the choice of ensemble for the MD calculations and which types of observables are best suited for each ensemble. In Sec. II B, we discuss the Kubo–Greenwood (KG) formalism for calculating the conductivity. In Sec. II C, additional optical properties beyond those requested in the TCCW2 are discussed. Optical properties of HED matter, such as the index of refraction and absorption, are relevant for experiments and can be calculated essentially for free once the conductivity is obtained. In support of this, we introduce a simple yet highly accurate algorithm for performing the principal value integration required by use of Kramers–Kronig relations. This is especially important when we look at the DC limit of various optical properties. Finally, we discuss transport coefficients that can be obtained from the Green–Kubo (GK) formalism, namely, the viscosity.

DFT-MD simulations typically operate in an ensemble with a fixed particle number (order 100–1000 atoms) with a fixed simulation cell volume. Invoking the Born–Oppenheimer approximation, the problem can be recast as ions moving in a potential determined by the electronic free energy corresponding to a given configuration of ions. Moreover, one typically assumes the ions can be treated as classical particles, which is almost always an excellent approximation in WDM regimes.

Calculating equilibrium properties in this ensemble is similarly straightforward. Through the use of a thermostat to control the temperature, we perform a molecular dynamics simulation to generate snapshots of ionic configurations, compute observables, and average the results over configurations. Energies, pressures, and electronic optical properties are examples of static observables amenable to this technique. Computing accurate dynamical properties within an NVT ensemble unfortunately requires much more care than directly using the NVE ensemble. The reason is because the method of enforcing the NVT ensemble in molecular dynamics calculations is through the use of a thermostat, which maintains a set temperature by imparting and removing energy from the ensemble of ions. This necessarily affects the ion trajectories and resulting time dependence of observables, though the impact can be more or less severe depending on the choice of thermostat.³¹ For the NVT ensemble, we use velocity rescaling to our target temperature as the thermostat.

As an example, in Fig. 1, we show the computed mean squared displacement (MSD) as a function of time for H in the priority 1 CH system. The MSD is directly related to the diffusion coefficients through the well known Einstein relation^32,33 and can also be obtained from the velocity autocorrelation function and is thus clearly a dynamical property. We find that both the short time behavior of the MSD and potentially even the slope (and therefore the self-diffusion constant) are noticeably different between the two ensembles.

We obtain accurate estimates of the dynamical properties in an NVT ensemble by drawing several statistically independent configurations from an NVT simulation, which are used as starting configurations for subsequent NVE simulations. We analyze the NVT simulation first to determine the autocorrelation length of properties such as the energy or pressure. Using that autocorrelation length, we sample atomic configurations that are spaced by at least that much. To ensure that we are fully equilibrated, the snapshots are obtained from the end of the NVT trajectory. We carry out NVE simulations on these snapshots and average their properties, which is consistent with the correct state sampling implied by Eq. (1). To determine if enough NVE samples are used, we ensure the average temperature and energies from NVE snapshots are consistent with the NVT values. Typically, 5–10 NVE snapshots are sufficient. Then, we use these same samples to calculate dynamic observables such as real time autocorrelation functions (and therefore quantities derived from those, such as viscosity, power spectra, diffusion coefficients, etc.).

When compared to experimental measurements, often one is most interested in the zero frequency, or DC conductivity, $σ DC = lim ω → 0 σ 1 ( ω )$⁠. In order to converge the electrical and thermal conductivity, there are a number of technical parameters that need to be addressed. While some of these technical details have been discussed previously in the literature,²⁷ we present a thorough discussion here. In particular, values need to be converged with respect to the number of Kohn–Sham states to include the number of k-points, the number of configurations from the MD trajectory, the finite size effect due to the size of the simulation cell, and the representation of the delta function.

In Fig. 2, we show the convergence of both the electrical and thermal conductivity as a function of the number of bands considered. A number of important features can be observed in Fig. 2(a). First, the increase in the number of states only contributes to the high frequency part of the conductivity. In many cases, only the DC conductivity is of interest, i.e., the $ω → 0$ limit, which rapidly converges and is not impacted by increasing the number of Kohn–Sham states. Note that the high frequency tail becomes important for other optical properties and thus may require significantly more bands. This is because the optical properties are functions of both σ₁ and σ₂, and σ₂ is related to σ₁ via an integral over all frequencies. The treatment of optical properties will be discussed further in Sec. II C.

In Fig. 3, we show the impact on the number of k-points. On a single snapshot, we calculate the KG expression by averaging over a variety of Monkhorst–Pack k-point grids. The dominant effect is to reduce the fluctuations and noise in the conductivity, which helps improve the estimation of the ω = 0 limit. We see that even moving from a 1 × 1 × 1 to a 2 × 2 × 2 Monkhorst–Pack grid, the magnitude of the fluctuations is significantly reduced. In practice, static DFT calculations are carried out on multiple uncorrelated configurations sampled from an equilibrated NVT trajectory. The conductivity is obtained from the average over configurations and is illustrated in Fig. 4. Averaging over multiple ionic configurations is required to produce a converged value of the electronic contribution to the conductivities, but the ionic contributions, expected to be much smaller, are here neglected. The DC conductivity rapidly converges with the number of snapshots, but the fluctuations and statistical uncertainty are reduced with more samples. Therefore, convergence in both the number of k-points and configurational averaging is needed to produce a converged estimate of conductivity.

There is an unphysical drop in the conductivity near zero frequency due to finite size effects. This comes from the discrete eigenvalue spacing due to including only finitely many ions.⁴¹ The artificial periodicity imposed on the ionic configurations by the choice of cell size introduces degeneracies in the eigenstates. Ionic density fluctuations $ρ k$ that are not commensurate with the simulation cell G vectors are not sampled, and so these degeneracies are not lifted, thereby introducing artificial gaps into the DOS and conductivity calculations. This results in a large finite-size effect in the DC limit of the conductivity. This can be seen visually by inspecting the eigenvalue spacing near the Fermi energy as these states are the dominant contributor to the conductivity.

We investigate the impact of using larger simulation cells for CH in Fig. 6. We see that for the smallest supercell, the conductivity in the 0.1–10 eV range is underconverged relative to the 256 and 500 atom simulation cells, which are in good agreement with one another. In addition, the low frequency behavior is extremely sensitive to the supercell size, which can be seen in both the electrical and thermal conductivities. By increasing the simulation cell size, the intercept at ω = 0 systematically increases; in the limit of an infinite simulation cell, we would expect to recover the correct DC limit. However, this convergence is quite slow in simulation cell size and would require even larger simulation cells in order to have a reliable estimate of the DC conductivity directly from the KG formalism.

There are some additional considerations for converging the KG conductivity. In order to evaluate the expression, an approximation is used to represent the delta function expression given in Eq. (3). Typical choices are Gaussians or Lorentzians; the Gaussian representation is the most common choice in the literature, whereas the Lorentzian is a physically meaningful approximation that reproduces the exact Kubo–Greenwood expression in the limit of zero smearing, as discussed in Calderín et al.⁴² The main impact of each smearing function and smearing width is to smooth out the fluctuations, as can be seen in Fig. 8. If the smearing parameter is too large, the incorrect DC limit could be obtained as the conductivity is overly smoothed, particularly using the Lorentzians. If the smearing is too small, the conductivity is significantly more noisy. Therefore, the smearing parameters are scanned to find a balance between reducing noise and guaranteeing the correct DC limit. For the data shown in Fig. 8, smearing values above 0.15 eV would yield biased estimates of the DC conductivity.

As discussed in Calderín et al.⁴² and described in Sec. II B, the standard expression for the Kubo–Greenwood conductivity does not include intraband contributions, which could be important for some materials. For $L m n α , β ( ω )$ in Eq. (3), we see that the eigenvalue difference in the denominator is divergent if the eigenvalues are the same. This can occur for both degenerate states as well as the contributions from the same state. Separating the sum into degenerate, intraband, and interband contributions, we can replace the occupation/energy ratio by a derivative of the Fermi–Dirac distribution for the degenerate and intraband contributions. We note that for all systems considered here, we have tested both with and without intraband contributions and saw no significant differences in the conductivity. Unless otherwise stated, values reported here include intraband contributions.

The real part of the electrical conductivity, $σ 1 ( ω )$⁠, is determined from the KG formalism as described previously. Due to causality, the imaginary part of the conductivity $σ 2 ( ω )$ can be obtained from the KK transformation given in Eq. (9). Since we only have a numerical representation of $σ 1 ( ω )$⁠, we have to numerically integrate this expression to obtain $σ 2 ( ω )$⁠. Many data analysis packages, such as SciPy, Mathematica, MatLab, etc. provide numerical implementations of the Cauchy principal value integral. However, we have implemented our own analysis software for KG in C++ and therefore need our own numerical implementation. Rather than incorporating a third party numerical integration library to handle this, we present our own approach, which rewrites the integral in a form that allows us to analytically integrate the divergent terms. We also compare the accuracy of this approach to other approaches on a Drude model system to demonstrate the accuracy.

We apply all of the methodologies and convergence analysis discussed in Sec. II to a subset of the systems requested in the TCCW2. For all calculations, we have used VASP v6.3.2.^24–26 For the matrix elements used in the KG expressions, we have modified VASP to report the longitudinal matrix elements in a format suitable to our post-processing tools. A brief summary of the systems, including the density, temperature, functional, pseudopotential, and cutoff, is shown in Table I. The materials at each of these conditions are in the liquid phase. For the single species systems, we start with a simple cubic conventional cell and multiply by a diagonal transformation matrix to obtain a cubic simulation cell with a target number of atoms. Then, we scaled the lattice parameter to obtain the desired density. For multi species systems, a similar method was used with alternating species on different lattice sites. The liquid structure is obtained by running the initial DFT-MD in the NVT ensemble with the target temperature and allowing the system to equilibrate.

We summarize our results in Table II. In addition to the quantities reported in the main TCCW2 paper⁵⁰ (P, η, κ, and $σ DC$⁠), additional results are presented for optical properties such as the index of refraction n, the extinction coefficient k, the absorption coefficient α, and the reflectivity R. For these properties, values are presented at wavelengths of 1550 and 532 nm, which are commonly used wavelengths for various HED experimental diagnostics. We note that while there are no experimental measurements to validate any of these predictions, the approach presented here results in values that are in excellent agreement with other ab initio predictions provided in the TCCW2 workshop.⁵⁰ We would like to point out that our results are also in good agreement with other DFT-MD results, which used finite-temperature functionals,^21–23 which are thought to be more accurate for materials in the WDM regime.

In this work, we have extensively discussed many aspects required for generating accurate and reproducible estimates of transport coefficients in the WDM regime using DFT-MD. While attending the TCCW2, it became clear that estimating transport coefficients is often an idiosyncratic process unique to individual researchers and groups. Therefore, we have focused on describing in detail how to systematically produce converged and accurate estimates of these quantities, an aspect of the work that has been poorly documented in the literature. We described how the choice of ensemble can impact the results, how different observables converge with respect to various adjustable parameters available in the DFT simulation, as well as present discussions on finite-size effects. This should hopefully remove potential sources of uncertainty in our data when trying to compare against other DFT-MD calculations, as well as make it easier for third parties to reproduce our data.

Finally, to support the accurate calculation of more general optical properties, we introduced a simple yet accurate method for performing principal value integration, which is required when using the Kramers–Kronig relations to calculate the full complex conductivity, dielectric function, indices of refraction, and reflectivity. We tested this method extensively and showed much greater accuracy in the low frequency regime than other commonly used methods.

This work was supported by the Laboratory Directed Research and Development program (Project No. 229428) at Sandia National Laboratories, a multimission laboratory managed and operated by the National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA0003525.

Sandia National Laboratories is a multi-mission laboratory managed and operated by the National Technology & Engineering Solutions of Sandia, LLC (NTESS), a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration (DOE/NNSA) under Contract No. DE-NA0003525. This written work is authored by an employee of NTESS. The employee, not NTESS, owns the right, title, and interest in and to the written work and is responsible for its contents. Any subjective views or opinions that might be expressed in the written work do not necessarily represent the views of the U.S. Government. The publisher acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this written work or allow others to do so, for U.S. Government purposes. The DOE will provide public access to results of federally sponsored research in accordance with the DOE Public Access Plan.

The authors have no conflicts to disclose.

Cody A. Melton: Data curation (equal); Methodology (equal); Software (equal); Writing – original draft (lead). Raymond C. Clay III: Data curation (equal); Methodology (equal); Writing – original draft (supporting). Kyle R. Cochrane: Data curation (equal); Methodology (equal). Amanda Dumi: Data curation (equal); Methodology (equal); Writing – original draft (supporting). Thomas A. Gardiner: Formal analysis (equal); Methodology (equal). Meghan K. Lentz: Data curation (equal); Methodology (equal). Joshua P. Townsend: Data curation (equal); Methodology (equal); Software (equal); Writing – original draft (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.