Yukawa Friedel-Tail pair potentials for warm dense matter applications

Studies in warm dense matter (WDM) systems straddle an intermediate regime of matter extending from cold solids to hot plasmas, at arbitrary compressions and states of ionization GaffneyHDP18 ; Grabowski20 ; Hungary16 ; cdw-cpp15 ; NgIJQC12 . The ions are usually modeled as classical particles whereas the electrons are treated quantum mechanically stanton2018multiscale . Typical WDM systems often fall into regimes of average mass density $\bar{\rho}$ and temperature $T$ where reliable experimental data do not exist or are not easily available, as in most astrophysical and geophysical applications. The same issue arises in high-energy-density applications where energy is deposited in femto-second time-scales producing systems containing multi-species quasi-equilibria associated with multiple temperatures NgIJQC12 ; CDW-cm-91 ; ELR98 . The electrons and holes in semi-conductor nanostructures also behave like WDM systems due to the small effective masses of the charge carriers lfc1-dw19 . The electronic and ionic interactions in these finite-temperature systems are such that interparticle effects are strong; hence simple models of the effective pair potential fail when applied to such many-particle systems Stanek21 .

Furthermore, the equations of state (EOS) GaffneyHDP18 and transport properties  Grabowski20 of WDM systems require the tabulation of many data points and the available many-body quantum techniques, e.g., density functional theory (DFT), path-integral Monte Carlo (PIMC) methods etc., become computationally expensive. Here we aim to leverage results of DFT based $N$-ion simulation codes VASP ; ABINIT to (i) validate the functional form of physically motivated models to represent pair interactions, and (ii) elucidate the regions of physical space for which simpler models are valid. Using validated pair potentials, we avoid carrying out repeated (“on-the-fly”) $N$-ion DFT calculations as the ion system is evolved using “classical molecular dynamics”. Here, “classical molecular dynamics” implies that the electronic structure calculation needed for ion interactions is pre-computed. We emphasize that this pair potential description does not neglect many-ion effects as explained below. The relationship between the pair energy and the pair potential, and how higher-order corrections are included, are considered in Sec. VII. The standard DFT-MD procedure that determines the electronic structure of an $N$-atom sample of the material via Kohn-Sham-Mermin finite-temperature DFT (e.g., as available in computational packages VASP ; ABINIT ), together with the use of classical MD to evolve the ions to equilibrium will be referred to here as “QMD” for brevity.

The computational cost associated with $N$-ion DFT can be reduced by using a single-ion DFT approach DWP82 ; eos95 ; Hungary16 ; cdw-pop21 where an ion-ion exchange-correlation (XC) functional is used to reduce the many-ion problem to that of a single ion such that the material is a simple superposition of pseudoatoms. This does not mean that a simple superposition of actual atoms is used, as may be possible in weakly coupled systems. This “single-ion” DFT approach to the many-ion, many-electron problem is referred to here as the neutral-pseudoatom (NPA) model. The NPA calculations, being similar to average-atom (AA) calculations, are computationally rapid in comparison to QMD calculations, and usually produce EOS data and pair potentials that are in good agreement with QMD results Stanek21 ; cdwSi20 ; Hungary16 . The NPA calculations also directly yield single-ion properties like the mean-value of the charge ($\bar{Z}$) carried by ions, electron-ion scattering cross sections, or electron-ion pseudopotentials that are useful for transport calculations Grabowski20 . Note that there are many AA models SternZbar07 ; Rosznyai08 ; FausAA10 ; Murillo13 ; StaSau13 ; StaSau16 that are similar in spirit but differ in their details. Examples include different choices of functionals used to treat exchange-correlation and kinetic energy, and more significantly, the boundary conditions used in defining the AA which may be confined to a Wigner-Seitz sphere or to a much larger volume DWP82 . The particular NPA model used here is described in detail in Refs. cdwSi20 ; cdw-Carbon10E6-21 ; eos95 and will be referred to here generically as “NPA.”

Recently, machine learning (ML) methods have been employed to accelerate the calculation of multi-center interatomic potentials (MCPs) whose simplest member is the pair potential. One of the simplest ML methods is force matching ercolessi1994interatomic which aims to minimize a loss function constructed from interionic forces from a QMD dataset and a predetermined force law Stanek21 ; Brommer2007Potfit ; Brommer2006Potfit ; Brommer2015Potfit . Such a loss function may take the form

where $M$ denotes the number of particle coordinate configurations, $F_{k}(\xi),E_{\ell}(\xi)$, and $S_{m}(\xi)$ are the force component, energy, and stress tensor component from the target interaction potential parametrized by $\xi$ with weights $w_{k}^{F},w_{\ell}^{E},$ and $w_{m}^{S}$. The reference force component, energy, and stress component, $F^{0}_{k},E^{0}_{\ell}$, and $S^{0}_{m}$, are determined from high-fidelity calculations such as Kohn-Sham DFT MD. The last two terms in Eq. (1) can be thought of as regularizers of the loss function that restrict the functional from of the optimized interaction potential. In the following, all of our force-matched (FM) potentials assume a loss function given by Eq. (1) with $w_{\ell}^{E}=w_{m}^{S}=0$. More advanced machine learning approaches allow for greater than pair or three-body interactions by including density-density correlations thompson2015spectral . Deep neural networks have also been employed to learn the QMD potential energy surface where translational, rotational, and permutational symmetries are enforced zhang2018deep ; DeepLearning-QZeng21 . In all of the aforementioned ML methods the MCP is obtained for use in classical MD simulations. In the MCP approach the electron sub-system is suppressed, and hence higher-order (e.g, 3-body, 4-body) potentials are needed to capture the electronic energies that are not explicitly included; the philosophy is similar to that of the classic Stillinger-Weber potentials SW85 .

In contrast, the NPA approach uses only pair potentials, but explicitly includes electron contributions to the total free energy directly available from the NPA calculation, as indicated in more detail in Sec. VII. The NPA uses the electron density to calculate energy contributions from the electron subsystem whereas the MCP method tries to embody them in the MCPs that carry “electronic bonding information.” This is claimed to allow the MCP method to deal with structural defects (e.g., dislocations in metals) “on the fly,” without resorting to explicitly handling electronic effects arising from such inhomogeneities. There is however no underlying principle in physics which asserts that there exists a general classical MCP that determines the physics of general electron-ions systems in the way claimed by ML; the only available principle, due to Hohenberg and Kohn, is that there exists a one-body density functional that is universal.

An advantage of ML methods is that the functional form of the ion-ion pair interaction potential is allowed to be determined by data, eliminating any bias associated with assuming too simple a functional form based on a limited physical model. On the other hand, the lack of a functional form prevents accurate recovery of asymptotic forms imposed by the physics of these systems. The data-driven nature of these ML methods is also one of their main limitations: the pair interaction potential will only be as accurate as the data used in the training process. It is possible that the QMD dataset may consist of too few time steps, resulting in a poor representation of the potential energy surface. The QMD dataset may also contain errors that originate from simulations of a small number of particles (referred to herein as finite-size effects). Quantifying errors from finite-size effects is critically important near phase transitions where the cooperative effect of a large number of particles occurs.

Furthermore, as very complex $N$-atom structures become possible with large-$N$ simulations, increasingly more advanced XC functionals become imperative. It has been shown that numerical results using different XC functionals at different levels of the “Jacob’s ladder” may give significantly different results in such calculations Remsing17 .

To mitigate the concerns associated with dataset inadequacies, we present a parametrized analytic form of a pair potential based on a general physical model that is found to be able to fit either to potentials extracted from QMD simulations or given by NPA calculations. If the QMD dataset is expected to be comprehensive enough to describe the potential energy surface, the analytic fit form will correct for finite-size effects at both small and large interatomic distances.

We note that the NPA is a DFT method where both the electron many-body problem and the ion many-body problem are jointly reduced to equivalent non-interacting one-body problems. This is done by relegating higher-order effects beyond mean-fields into exchange-correlation (XC) functionals; these depend only on the electron one-body density and the ion one-body density. Since a central ion is used as the origin of the NPA calculation, the one-body ion-density becomes effectively a two-body ion density and correlations beyond the two-body potential appear in the ion-ion XC-correlation functionals; classical ions have no exchange, and the correlations used in the NPA for solids, liquids, or hot plasmas are those brought in by the Ornstein-Zernike procedure. Thus, while the ion-ion XC-functional is strongly non-local, the electron-electron XC-functional used in the NPA is the finite-$T$ functional given by Perrot and Dharma-wardana PDWXC , within the local-density approximation. Due to the simplicity of the electron density profile in the NPA (which is a single-center screened ion), the local-density approximation is found to be sufficient.

In the following we treat simple fluids like liquid Al as well as complex liquids like liquid C, not only at high temperatures, but also at low temperatures where the static structure factors $S(k)$ and pair-distribution functions (PDFs) $g(r)$ may be thought to depend on bond-angle distributions and other complex correlation effects. However, density functional theory is a theory of the ground state energy for systems at $T=0$ and a theory of the free energy for finite-$T$ systems. These depend only on the time-averaged picture of the fluid which is uniform as we are not treating nanostructures or solids. Hence there is no need to include the detailed instantaneous bond-orientational details that are generated in standard $N$-ion simulations, only to be averaged out using millions of such calculations. The static (i.e., long-time average) correlations beyond the two-body terms are brought into the NPA calculations via the Ornstein-Zernike procedure which requires the use of an ion-ion XC functional, as shown in Ref. DWP82 . The energy effect of such correlations is taken care of in the NPA through the ion-ion XC functionals, as discussed in Ref. DWP82 . Furthermore, the energy and free energy are functionals of just the static one-body densities. Unlike in $N$-ion DFT simulations, calculations using the NPA method will not yield instantaneous bonding details and “snap shots” of the fluid, although the static PDFs and static structure factors become readily available. However, as shown by Hafner Hafner89 for liquid As, a simple pair potential can be used to provide considerable information on bond-angle distributions by using the pair potential in a classical simulation.

In this study, we focus on the pair potentials obtained from the NPA model, and on force-matched potentials extracted from QMD where they are available. They are of interest for simulating physical systems such as liquid metals, WDM systems, and strongly coupled plasmas where the primary input to classical MD simulations is the pair potential cdw-aers83 ; chenlai92 ; Utah12 ; VorbGeri-Pots13 ; whitley15 ; Stanek21 . Specifically, we use the NPA model to generate pseudopotentials and pair potentials needed for computations of physical properties of finite-temperature electron-ion systems of uniform density, and show that a simple physically meaningful parametrization of the potentials exists for them. The simplest method, based on a linear-response approach to the pseudopotentials, is found to be applicable for a wide range of densities and temperatures. Examples where linear-response fails are common for expanded metals, liquid (which we denote as $l$-) C, and $l$-transition metals at low $\bar{\rho}$ and $T$ Stanek21 . In such cases special procedures and sophisticated XC-functionals are needed, even in standard $N$-atom DFT.

In this study, we provide pair potentials for simple fluids like $l$-Al and complex fluids like $l$-C whose QMD simulations show the existence of transient covalent bonds. Nevertheless, it was shown in detail that a simple-metal model is applicable to the thermodynamic and structural properties of $l$-C and $l$-Si, and even to their phase transitions cdw-carb22 ; cdwSi20 . These results may seem to go against “chemical intuition” but are consistent with well established physics-based insights on the ionic correlations in simple and complex liquid metals and solids Hafner87 .

The remainder of this manuscript is organized as follows. First, in Sec. II, we present a brief overview of the NPA model and pair potentials. We discuss the theory of NPA pseudopotentials which originates from the Kohn-Sham DFT electron density. We also describe how pair potentials stem from linear-response theory and their implementation in the NPA model. The NPA pair potentials are shown to agree closely with force-matched QMD pair potentials which are available for length scales limited to the size of the simulation box used in QMD. Then, in Sec. III, we introduce a simple but generally applicable parametrization scheme for pair potentials in homogeneous metallic systems. The proposed parametrization is consistent with a simple physical picture based on screened Coulomb interactions of the Yukawa type, and includes thermally damped Friedel oscillations as they become long-ranged and important for $T<E_{F}$. In Secs. IV to VI we present fits to our parametrization scheme using pair-potentials from NPA calculations for $l$-C, $l$-Al, and $l$-Li. We show that our parametrized scheme accurately reproduces structural quantities obtained from QMD or PIMC calculations and can be used to greatly reduce errors incurred from finite-size effects in QMD simulations.

The NPA calculation self-consistently generates the Kohn-Sham one-body electron density distribution $n(r)$ around a nucleus placed in the appropriate environment in the fluid where an all-electron calculation is used. For simplicity of exposition here we consider the case where the bound electrons are well localized inside the Wigner-Seitz sphere of radius $r_{\rm ws}$ associated with the ion. Then, the electron distribution can be unambiguously written as a sum of a bound electron part denoted as $n_{b}(r)$ and a free electron part denoted as $n_{f}(r)$, while the issues arising from the more general case are discussed in Appendix A. The free electrons are not confined to the Wigner-Seitz sphere (as in many AA models), and occupy the whole of space represented by a large “correlation sphere” of radius $R_{c}$ which is some five to ten times $r_{\rm ws}$, chosen to ensure that all interparticle correlations are sufficiently small by $r=R_{c}$, i.e, all pair distribution functions (PDFs) $g(r)\to 1$ as $r\to R_{c}$. The electron chemical potential is the non-interacting value used in DFT and need not be re-assigned as in many AA models that use a Wigner-Seitz sphere to define and AA. The free-electron density $n_{f}(r)$ is used to construct a pseudopotential which is then employed to construct a pair potential. If a metallic model with strong screening is assumed, the above steps can be carried out within linear-response theory, giving rise to simple local (i.e., $s$-wave) pseudopotentials that contrast with the non-local non-linear pseudopotentials supplied with QMD codes VASP ; ABINIT .

At the end of an NPA calculation we have the free-electron density distribution $n_{f}(r)$ around the ion of mean charge $\bar{Z}$, where the average (uniform) free-electron density $\bar{n}=\bar{Z}\bar{\rho}$ as required by charge neutrality. The evaluation of $\bar{Z}$ used in the NPA is both an atomic physics calculation and a thermodynamic ionization-balance calculation that satisfies the Friedel sum rule and charge neutrality.

The ion-ion pair potential, given in Eq. (5), when Fourier transformed gives a numerical form for $V(r)$. Since various model forms based on the Yukawa potential and its modifications have been studied in the context of WDM, it is of interest to examine the structure of $V(r)$ and provide a general physical model that can successfully parametrize these potentials accurately, using only a minimum of parameters that also have a physically clear meaning.

Such a model consists of a Yukawa-screened ion carrying a long-ranged Friedel tail referred to in the following as the “YFT” model. Here, we provide parametrizations for Al, C, and Li as examples of YFT potentials. We find that the model is useful even in regimes where the NPA technique becomes inaccurate, for example, for low-density low-temperature $l$-C, as shown by Stanek et al. Stanek21 . The FM low-density low-temperature potentials obtained from QMD simulations using Eq. (1) do not recover many Friedel oscillations due to finite-size effects. However, the FM potentials can be reliably extended by fitting to the YFT model explicitly given in Eq. (8).

The YFT form of the potential used is taken as a sum of a Yukawa form $V_{\rm y}(r)$ and a thermally damped Friedel tail $V_{\rm ft}(r)$

Here the Yukawa form has a screening constant $k_{\rm y}$ which is of the order of the finite-$T$ Thomas-Fermi screening length, while a different $T$ and $\bar{n}$ dependent screening constant $k_{\rm ft}$ defines the thermally damped Friedel oscillations. The Friedel oscillations have the usual form: a circular function phase-shifted by $\phi_{\rm ft}$ and modulated by a $1/r^{3}$ decay. This six-parameter fit form is found to be capable of accurately reproducing the pair potentials of all the WDM systems containing free electrons and ions that we have studied ranging from very low temperatures close to the melting point, to temperatures typical of the plasma state. At higher temperatures and higher compressions, the interaction potentials become simpler and Yukawa-like, with the Friedel oscillations becoming increasingly damped Stanek21 . Hence the potentials discussed here focus on the more demanding regime of WDM systems, i.e., from close to the melting temperature to about $\theta\sim 1$. In the Yukawa-like high-$T$ high-density regimes, Thomas-Fermi models and their extensions like the orbital-free molecular dynamics method can increasingly supplant $N$-atom DFT calculations. However, even in the high-$T$ range, the NPA provides a wealth of information (via its Kohn-Sham states), ion-ion PDFs etc., that are not available from Thomas-Fermi models or their extensions to “orbital-free” techniques, while being computationally very efficient.

The six parameters $a_{\rm y},k_{\rm y},a_{\rm ft},k_{\rm ft},q_{\rm ft},$ and $\phi_{\rm ft}$ can only be approximately estimated using analytical models of screening. This is partly because they are also playing the role of fitting to the form factors contained in the actual NPA potentials that account for the finite size of the ionic cores. Hence they are best obtained by fitting the NPA-derived numerical pair potentials to the YFT fit form Eq.(8). The simple least-square fitting routines available in standard open-source software or graphics packages (e.g., gnuplot or python) are sufficient for obtaining the fit parameters of Eq. (8). The fitted forms can then be used as inputs to classical MD simulations. For example, these fitted potentials can be used in systems with density gradients that can be described by a graded stack of layers of uniform density and width not less than the radius of the correlation sphere, viz., $R_{c}\sim 5r_{\rm ws}$.

Furthermore, the fitted pair potentials can be compared with the pair potentials derived from force matching to QMD simulations if FM potentials are available. Since $N$-atom QMD calculations are carried out in a simulation box of linear dimension $\sim r_{\rm ws}N^{1/3}$, the $N$-center potential energy surface used in force matching is itself of restricted range, and long-ranged oscillations are only poorly recovered in force-matched potentials unless a large number of particles are used. We illustrate this by comparing the QMD-derived C pair potential shown in Fig. 3(a) due to Whitley et al. whitley15 and the Li pair potential shown in Fig. 4 due to Stanek et al. Stanek21 . Such QMD-derived pair potentials can be approximately extended to include their long-ranged Friedel oscillations in two ways: ( i) they can be augmented with the Friedel tail obtained from YFT fits to NPA pair potentials or (ii) the YFT form can be fit directly to the QMD data. The parametrized forms can be analytically differentiated or manipulated more efficiently and more accurately than numerical tabulations or ML potentials.

Liquid C is usually considered to be a highly covalently bonded system even in the fluid state. This view is supported by “snap-shots” of real-space positions of the ions obtained from QMD simulations. The multi-center bond-order potentials Abell85 ; Brenner02 ; ghiring05 use this chemical picture to model MCPs for C systems. However, as is evident from DFT, only the long-time average is relevant to determining the thermodynamics and static structure factor since the effects of the transient bonds average out in $l$-C DWP-carb90 . In QMD, this average is explicitly carried out using millions of configurations of solid-like ionic clusters, without taking advantage of the spherical symmetry of the fluid or the possibility of using an ion-ion XC-functional to treat many-center effects. While the chemical-bonding picture does not use the Fermi-liquid character of $l$-C, the NPA approach exploits it to the fullest extent in the YFT model.

We calculate the NPA pair potentials and the $g(r)$ of $l$-C, at density 12.64 g/cm³, at $T=8.62$ eV (1x$10^{5}$K), 21.5 eV (2.5x10⁵K), and 430.886 eV (5x10⁶K) to compare with the PIMC and QMD calculations of Driver and Militzer driver12 . The mean ionization $\bar{Z}$ obtained from the NPA changes from $\bar{Z}=4.0$ at low $T$, to $\bar{Z}=4.05$ by $T=100$ eV, and reaches $\bar{Z}=5.758$ by $T=430.9$ eV – the highest $T$ used by Driver and Militzer driver12 .

The YFT fits to the NPA potentials that reproduce the PIMC $g(r)$ data even at low-$T$ and at lower density are given in Fig. 5. The comparison of the NPA $g(r)$ obtained from the NPA potentials with the $g(r)$ from PIMC are given in Fig. 6. The NPA $g(r)$ are obtained from the pair potentials using the MHNC equation. This is an inexpensive alternative to carrying out an MD simulation. We caution that if using MHNC in place of MD, one needs to specify an accurate bridge function in the MHNC equations. Additionally, although hard-sphere approximations to the bridge function are often employed in the MHNC equations by appealing to a universality principle Rosenfeld93 , it should be used with some care as this model may not be valid for certain systems. For example, the bridge function can be set to zero for C, Si, and Ge in the liquid state at normal densities, and near their melting points when their structures are controlled mainly by the scattering of electrons by ions from one edge of the Fermi surface to the opposite edge, giving rise to $2k_{F}$ scattering effects cdw-carb22 that control the ionic structure, rather than packing effects.

In Table 1, we present a sample of parametrizations for the YFT potentials for C along the $T=1$ and 2 eV isotherms, for the range of densities $\bar{\rho}=2.267$ g/cm³ up to 10 g/cm³. At higher $\theta$ and higher densities, the potentials increasingly approximate a Yukawa-like form. The linear response potentials are applicable for $\bar{\rho}\geq$ 3 g/cm³ where they recover the structure factor and pressure obtained via QMD simulations cdw-carb22 . The low-density range from $\bar{\rho}<3$ g/cm³ and extending towards the graphite density and below, cannot be accurately computed using the linear-response model given by Eq. (2). Therefore, we use a FM pair potential to construct a parameterization with the YFT model and provide parameters for the density $\bar{\rho}=2.267$ g/cm³. We note that the mean ionization of C is such that $\bar{Z}=4$ for the densities and temperatures studied here.

Higher densities and temperatures (e.g., 10-100 times the diamond density and in the 10⁶ K temperature range) have been studied using the NPA cdw-Carbon10E6-21 , and the pair potentials under such conditions can be accurately fitted to the YFT form. To illustrate the fact that the YFT potentials reduce to simple Yukawa-like potentials for sufficiently high-$T$, we give the parametrizations for the Yukawa-like pair potentials of $l$-C at the graphite density in Table 2 of Appendix B.

The parameters given in Table 1 reflect the character of a screened metallic ion. This is confirmed by the value of the screening constant $k_{\rm y}$ which is quite close to the Thomas-Fermi screening wavevector, while the wavevector $q_{\rm ft}$ associated with the Friedel oscillations, is approximately 2$k_{F}$; the Friedel oscillations decay as $1/r^{3}$ as expected for a metallic fluid. The screening constant $k_{\rm ft}$ introduces a damping of the Friedel oscillations. For systems such that $\theta$ is small, this damping constant may be set to zero and the Friedel oscillations are similar to those familiar in metal physics. Since the $l$-C systems studied here are such that $\theta$ is very small, we could have set $k_{\rm ft}=0$ and used a simpler five-parameter fit which is equally capable of fitting to the QMD or NPA data. We will make this simplification for treating $l$-Al as is done in Sec. V and presented in Table 3.

We study normal density Al from $T=0.5$ eV to 100 eV, and in the density range $\bar{\rho}=$ 2.0 g/cm³ to 6 g/cm³ at $T=1$ eV. The normal solid-state density $\bar{\rho}$ is 2.7 g/cm³, while the isobaric $\bar{\rho}$ is 2.375 g/cm³ at its melting point ($T=$ 933.47 K). Unlike the C ion with $\bar{Z}=4$ which is almost point-like – having only the 1$s$ of bound electrons – Al³⁺ has a large core and its form factor is an important feature of its pseudopotential. Nevertheless, we find that the YFT model works satisfactorily for $l$-Al as well because the YFT parameter set is flexible enough for doing “double duty”. What we mean by this is that the parameters of the YFT model have more fitting capacity than would be found in a strictly defined physical model. Such a strictly defined model would require a parametrization of the form factor of the pseudopotential and a parametrization of the electron response function, so that Eq. 5 is itself fully parametrized.

A highly non-local pseudopotential rather than the simple $s$-wave form has been used in linear-response DRT75 and in QMD codes like the VASP VASP and ABINIT ABINIT although we find this to be unnecessary in NPA applications to uniform-density systems and cubic solids. Liquid Al differs from $l$-C in that the main peak of the structure factor is not located at 2$k_{F}$ as in C, Si, and Ge cdw-carb22 at their typical densities. The YFT parameters that fit the NPA potentials for $l-$Al are given in Table 3 of Appendix B. The WDM form of $l$-Al very close to its melting point, and under more extreme conditions have been studied using the NPA and reported in previous publications eos95 ; cond3-17 ; PDWBenage02 .

The Friedel oscillations in $V(r)$ are of great importance for the potential near the melting point of Al. Under these conditions, the first peak of the $g(r)$ locates at a positive energy turning point of a Friedel oscillation DRT75 ; cdw-aers83 ; HarbEOSPhn17 ; DSF18 . The NPA potentials between $T=$0.5 eV and 9.0 eV have been compared with FM potentials obtained using Eq. (1) from QMD simulations Stanek21 . In Fig. 2 we show two cases of Al at $\bar{\rho}=2.7$ gcm^-3, viz., one for $T=$ 0.5 eV and another for $T=$ 9.0 eV. The $T=$ 0.5 eV case is typical of the YFT regime where Friedel oscillations are non-negligible, whereas the $T=$ 9.0 eV case typifies the regime where the Friedel oscillations have become sufficiently damped and the YFT potential is approximately Yukawa-like.

At low-$T$, the NPA potentials have oscillatory tails while the FM potentials capture them only partly, as already seen in the C potentials in Fig. 3(a). However, the FM pair potentials and the more oscillatory NPA potentials nevertheless yield PDFs and ion-ion structure factors that are in close agreement. By $T=9.0$ eV, which corresponds to $\theta=0.77$, the potentials are no longer oscillatory and only the Yukawa form remains. This is in fact often the case by $\theta\gtrsim 2/3$ even at other densities, and for many elements. As $T$ increases, $\theta$ increases and the Yukawa approximation improves, but this is modified by changes in the mean ionization $\bar{Z}$ which increases $E_{F}$, thereby countering the accuracy of the Yukawa form. Nevertheless $\theta$ steadily increases in WDM Al from $\theta=1.549$ for T = $20$ eV where $\bar{Z}$ = 3.495 to $\theta=4.567$, for T = 100 eV where $\bar{Z}$=7.721. Thus the Yukawa form continues to be a reasonable approximation in spite of increased ionization; the values of the parameters for the range $T=8$ eV to 100 eV for normal density Al is given in Table 4 of Appendix B.

Although the potentials have become Yukawa-like, the ion-ion $g(r)$ and $S(k)$ are still strongly coupled as the plasma parameter $\Gamma=\bar{Z}^{2}/(r_{\rm ws}T)$ exceeds unity, being $\Gamma=8.35$ at $T=10$ eV and $\Gamma=6$ at $T=80$ eV. Hence $g(r)$ must be computed using MHNC or with MD. A simple but useful approximation may be obtained by a one-component-plasma analogy based on the value of $\Gamma$ Clerouin14Gamma . The fact that the potentials reduce to Yukawa-like forms in at least the $\theta\gtrsim 1$ regime opens the door for a semi-analytic theory of these WDM systems for a wide range of conditions.

In fitting the NPA potentials to the YFT form we have simply used no thermal damping on the Friedel term ,i.e., setting $k_{\rm ft}$ equal to zero. This shows that the parameters of the fit form “do double duty” and their values are not necessarily related to the values expected from basic theoretical models, although quite close to them. This fortunate flexibility is partly associated with the freedom in the selection of the range of $r$ used in the fits. Standard pseudopotential formulations also have such flexibility arising from how the atomic core is excluded in the pseudopotential. In general, the fit range is selected to be that of the ion-ion $g(r)$, viz $r_{\rm ws}\leq r<R_{c}$.

It should be noted that the pair potential defined by Eq. (5) does not contain interactions with the polarizable atomic cores of the two interacting Al ions, or van der Waals interactions. The YFT fit-forms also do not specifically contain such interacting terms, although the parameters have some flexibility to absorb such contributions approximately. How these concerns could be included within a purely DFT scheme were discussed in the Appendix of Ref. eos95 . These concerns are also important in systems that contain a significant mixture of neutrals, e.g., Ar plasmas at $T\sim$ 1$-$2 eV and normal density Stanek21 where an argon mixture is treated using g the NPA. At the temperatures and densities considered here, the dominant interactions are those included in Eq. (5), and hence we will not discuss core-polarization effects further.

In Fig. 9(a) we compare the $g(r)$ generated directly within QMD, and labeled QMD-VASP with the $g(r)$ obtained from MD using the FM potential extracted from the QMD using Eq. (1). In Fig. 9(b) we compare the $g(r)$ generated directly from QMD using ABINIT, with the PDFs obtained from the NPA pair-potential using MD, and using MHNC DSF18 , showing that the MHNC yields $g(r)$s of good accuracy. We may also note that the NPA pair potential provides an $S(k)$ for Al that not only agrees with the static $S(k)$ obtained by QMD simulations, but also accurately renders the dynamic ion-ion structure factor $S(k,\omega)$ obtained from DFT-MD microscopic simulations available at $T=$ 1 eV and at $\bar{\rho}$ = 2.7 g/cm³ DSF18 .

Lithium in its WDM state has not been as extensively studied as Al or C. However, accurate data for the structure factor $S(k)$ of $l$-Li near its melting point as well as attempts at theoretical predication are available Salmon2004 . Since theories should preferably be tested against actual experimental data rather than against simulation results, we first study this case, i.e., $l$-Li at $\bar{\rho}$ = 0.513(4) g/cm³ and at a temperature of $T=$ 197^∘C (which corresponds to $T=0.0405$ eV in the units used here). At this very low temperature and nearly the solid-state density, it is easy to imagine that pair potentials would be totally inadequate and that a multi-center approach would be needed. However, our results show that the NPA pair potential, taken together with the MHNC integral equation where the bridge term is modeled by a hard-sphere fluid very successfully predicts the experimental data. The packing fraction $\eta=0.3$ of the hard-sphere fluid is determined using the Lado-Foiles-Ashcroft free-energy minimization criterion LFA83 .

Stanek et al. Stanek21 have compared the NPA potential and its predictions for Li at $\bar{\rho}=0.513$ g/cm³ and $T=0.054$ eV, with those from FM potentials obtained from QMD simulations using Eq. (1), obtaining excellent agreement. As the Fermi energy of Li at this density is $E_{F}\sim$ 5 eV, with $\theta\simeq 0.011$, the electron subsystem is effectively highly degenerate. The NPA potentials were found to work very well even for solid-density Li as the phonons calculated using the NPA pair potentials were in good agreement with experimental phonon dispersion data HarbEOSPhn17 . It has also been verified that the NPA pair potential for Li at $T=0.173$ eV generates PDFs in good agreement with the QMD data at $\bar{\rho}=0.85$ g/cm³ Kietzman08 .

While the above cases are for high-electron degeneracy, Harbour et al. cond3-17 found good agreement between QMD calculations of various aspects of X-ray Thomson scattering and those from NPA methods for the density $\bar{\rho}=0.6$ g/cm³ and $T=4.5$ eV. In this case, the electron system is at a temperature close to the Fermi energy and is partially classical in its behaviour. Hence we may conclude that currently available QMD data confirm the validity of the NPA potentials for WDM-Li from the solid state to the plasma state. A comparison of the FM Li-Li pair potential using Eq. (1) at $T=0.54$ eV and $\bar{\rho}=0.513$ g/cm³ and the corresponding NPA Li-Li potential is given in Fig. 11.

Unlike in the case of Al, the difference between the FM $V(r)$ and the NPA $V(r)$ at their minima is more significant compared to the case of Al at $T=$ 1 eV (Fig. 2), being of the order of the thermal energy of the system at $T=0.054$ eV. However, this difference is not reflected in the predicted properties, e.g., the $S(k)$ or $g(r)$ Li, as seen in Fig. 12 and Fig. 13 respectively. In Fig. 12 we show the $S(k)$ obtained from the HNC, i.e., MHNC with the bridge function set to zero as well as the MHNC with the proper LFA-optimized bridge term LFA83 , displaying the role of the bridge contribution to the ionic structure.

In Table 5 of Appendix B we provide YFT fits to Li-Li pair potentials generated using the NPA electron density and the linear-response pair potential given by Eq. (5). In most cases, we have set the parameter $k_{\rm ft}$ controlling exponential damping term of the Friedel oscillations to zero, except in the case of the density $\bar{\rho}=$ 0.85 g/cm³ where we have included it as a parameter, merely to illustrate the fact that there is significant interplay among parameters, and consequently, a numerical fit of similar quality can usually be obtained even if $k_{\rm ft}$ is set to zero. We give YFT fits for several densities from $\bar{\rho}=0.513$ g/cm³ to 1 g/cm³, and at temperatures $T=0.054$ eV and 1.0 eV in all cases. We also give the YFT fit at $T=4.5$ eV and $\bar{\rho}=$ 0.513 g/cm³ which corresponds to $\theta=0.983$ where the fit form reduces to a Yukawa form as the Friedel oscillatory term is found to be negligible. Just as in Al, we have a regime $\theta\gtrsim 1$ where the pair-potentials of the Li WDM system become Yukawa-like in spite of increasing ionization with increasing temperature. The Yukawa-like parameters that fit the NPA potentials in this regime are given in Table 6 of the Appendix B.

There are multiple motivations for constructing parametrizations of pair potentials based on an approximate physical model with very few parameters rather than using a highly flexible purely numerical approach fitted to QMD data. First, the parametrization provides some validation to the associated physical picture. Second, the fits can be used in simulation codes which call for repeated evaluations of the potentials and total energies as a function of position, temperature, and density. Third, the model potentials can be used to extend the domain of force-matched potentials $V(r)$ available only for values of $r$ restricted by the size of the simulation cell. Fourth, the existence of large regions of $\bar{\rho},T$ phase space where the Yukawa form of the pair potential raises the possibility of a semi-analytic theory of a wide range of WDM systems that are treated by AA, orbital-free-Thomas-Fermi type calculations.

The existence of a simple YFT-type parametrization of pair potentials in these WDM systems supports the physical picture of them being essentially metallic liquids; this being true even in the case of $l$-C or $l$-Si where a different paradigm has often been advanced, within the intuitive chemical picture of covalent bonding. The paradigm of “bond-order potentials” and their limitations in WDM applications are now well known. Just as the “bond-order potentials” may be used in simulation codes that call for repeated evaluations of the potentials, these YFT potentials may also be similarly used, much more simply. The method based on constructing FM potentials extracted from more fundamental DFT-MD simulations may become computationally prohibitive especially at higher temperatures, although the underlying physics for $\theta\gtrsim 1$ is Yukawa-like.

Standard DFT fails when ground states contain significant non-Born-Oppenheimer corrections, while excited states may involve conical intersections. Nevertheless DFT-MD is an excellent method for testing methods like the NPA, or in providing inputs to reference data bases for systems where the Born-Oppenheimer approximation holds and when experimental data are not available.

We have presented an analytic form for ion-ion pair potentials for cold-, warm- to hot-dense-matter systems having a significant degree of ionized electrons, as found in metals and plasmas. Our model uses a Yukawa-like pair potential augmented by a thermally damped oscillatory Friedel tail, where the parameters in the potential are obtained from a single-center Kohn-Sham calculation within the NPA model or from QMD calculations. The YFT model is applicable at arbitrary coupling, unlike the standard Yukawa form valid only in the high-temperature high-density limit such that $\theta\gtrsim 1$. We have shown that it accurately reproduces results from QMD simulations for a wide range of plasma conditions from near the melting point to high temperatures at densities typically found in liquid metals and plasmas of interest in many areas of study.

Using reliable pair interaction potentials such as those from force matching to $N$-atom DFT-MD data, we find that in all cases studied here, the YFT agrees closely with them and reproduces structural quantities from DFT-MD with good accuracy. The YFT model is based on a linear pseudopotential that is consistent with the NPA calculation as well as its sophisticated ionization-balance calculations that go into the evaluation of the mean ionization per atom. The latter is given by $\bar{Z}$ and is self-consistently constrained to satisfy the Friedel sum rule. The pseudopotential and the mean ionization $\bar{Z}$ are primary inputs to transport calculations that are not easily available from $N$-atom DFT-MD or from orbital-free DFT calculations. The computational costs of NPA calculations, or YFT fits are orders of magnitude less than for QMD where millions of $N$-center DFT calculations for millions of ionic configurations are needed. We also note that in contrast to $N$-atom orbital-free molecular dynamics, NPA calculations (i) are computationally faster, (ii) do not suffer from finite-size effects, (iii) decrease statistical errors associated with short-time simulations, (iv) statistically resolve all species in systems with impurities, (v) are computable at any temperature, and (vi) generate pseudopotentials from their own all-electron calculations instead of requiring them as external inputs.

In the context of FM potentials, the YFT model provides a way to correct for finite-size effects; since it can be evaluated at small and large interparticle distances where as extrapolation is needed for FM potentials obtained from finite-cell calculations.

While the NPA model in its current formulation assumes the existence of a weak pseudopotential for which linear response theory is applicable, the YFT model is able to also accurately represent potentials obtained by force matching to systems like low-density $l$-C where NPA methods prove inadequate. The NPA fails if the system has a low density of free electrons (e.g., expanded metals and attenuated plasmas) or when no free-electrons what so ever (e.g., as in liquid molecular hydrogen) are present. In such situation, physics-informed ML methods are useful as they generate a very detailed description of the physical system including ion-ion many-body interactions and bonding descriptions inclusive of angular correlations. However, results from such ML methods as well as $N$-atom DFT methods should be validated against experimental data or strictly ab initio methods that directly obtain the wave function, as in configuration-interaction calculations, coupled-cluster calculations or quantum Monte Carlo and path-integral methods, when they are feasible.

While this work is focused on single-component systems, a natural extension of this work is to study plasma mixtures of different elements, or mixtures of different charge states of the same element eos95 . The possibility of treating mixtures would allow for the calculation of the equations of state, charge states of ions for plasma spectroscopy, and plasma transport properties like interdiffusion coefficients, that are of interest in hydrodynamic models used in simulations of high-energy-density experiments.

The authors would like to thank Raymond C. Clay III for useful discussions. L.J.S. and M.S.M. acknowledge support from the U.S. Air Force Office of Scientific Research Grant No. FA9550-17-1-0394.