Macroscopic simulations of dense plasmas rely on detailed microscopic information that can be computationally expensive and is difficult to verify experimentally. In this work, we delineate the accuracy boundary between microscale simulation methods by comparing Kohn–Sham density functional theory molecular dynamics (KS-MD) and radial pair potential molecular dynamics (RPP-MD) for a range of elements, temperature, and density. By extracting the optimal RPP from KS-MD data using force matching, we constrain its functional form and dismiss classes of potentials that assume a constant power law for small interparticle distances. Our results show excellent agreement between RPP-MD and KS-MD for multiple metrics of accuracy at temperatures of only a few electron volts. The use of RPPs offers orders of magnitude decrease in computational cost and indicates that three-body potentials are not required beyond temperatures of a few eV. Due to its efficiency, the validated RPP-MD provides an avenue for reducing errors due to finite-size effects that can be on the order of $\u223c20%$.

## I. INTRODUCTION

High energy-density science relies heavily on computational models to provide information not accessible experimentally due to the high pressure and transient environments. The plasmas in these experiments typically contain strongly coupled ions and partially degenerate electrons, which constrains our microscopic modeling choices to molecular dynamics (MD) and Monte Carlo approaches. Interfacial mixing in warm dense matter, for example, requires costly, large-scale MD simulations; but such simulations reveal previously unknown transport mechanisms.^{1} It is therefore crucial to quantify the efficacy of computational models in different regions of species-temperature-density space so that the cheapest accurate model can be exploited to address such problems.^{2–4} While it is desirable to use short-range, radial, pair potentials (RPPs) to maximize the length and time scales, *N*-body energies may be required in some cases. Few studies have been carried out that comprehensively assess the limitations of RPPs and the regimes of utility for the extant forms; given that the force law is the primary input into MD models, it essential to have quantitative information about these force models.

A wide variety of RPPs have been developed for modeling dense plasmas. In some cases the accuracy of the model can be inferred from its theoretical underpinnings; in other cases, comparison to higher-fidelity approaches or experiments is needed. Limitations of the RPP approximation are generally unknown unless compared to an *N*-body potential simulation result such as Kohn–Sham density functional theory molecular dynamics (KS-MD). Both KS-MD simulations and this comparison are time-consuming processes that are limited to the temperature regime in which the pseudopotentials necessary for KS-MD are valid.^{5,6} Moreover, comparisons between RPP-MD and KS-MD are limited in the literature, have not been carried out for a range of elements and temperatures, and are often validated with integrated quantities where individual particle dynamics have been averaged and results are subject to cancelation of errors.

In this work, we carry out KS-MD simulations for a range of elements, temperatures, and densities, allowing for a systematic comparison of three RPP models. While multiple RPP models can be selected,^{7–11} we choose to compare the widely used Yukawa potential, which accounts for screening by linearly perturbing around a uniform density in the long-wavelength (Thomas–Fermi) limit, a potential constructed from a neutral pseudo-atom (NPA) approach,^{12–15} and the optimal force-matched RPP that is constructed directly from KS-MD simulation data.

Each of the models we chose impacts our physics understanding and has clear computational consequences. For example, success of the Yukawa model reveals the insensitivity to choices in the pseudopotential and screening function and allows for the largest-scale simulations. Large improvements are expected from the NPA model, which makes many fewer assumptions with a modest cost of pre-computing and tabulating forces. (See the Appendix for more details on the NPA model.) The force-matched RPP requires KS-MD data and is therefore the most expensive to produce, but it reveals the limitations of RPPs themselves since they are by definition the optimal RPP.

Using multiple metrics of comparison between RPP-MD and KS-MD including the relative force error, ion–ion equilibrium radial distribution function *g*(*r*), Einstein frequency, power spectrum, and the self-diffusion transport coefficient, the accuracy of each RPP model is analyzed. By simulating disparate elements, namely, an alkali metal, multiple transition metals, a halogen, a nonmetal, and a noble gas, we see that force-matched RPPs are valid for simulating dense plasmas at temperatures above fractions of an eV and beyond. We find that for all cases except for low temperature carbon, force-matched RPPs accurately describe the results obtained from KS-MD to within a few percent. By contrast, the Yukawa model appears to systematically fail at describing results from KS-MD at low temperatures for the conditions studied here validating the need for alternate models such as force-matching and NPA approaches at these conditions.

In Sec. II, we discuss how RPPs arise from second order perturbation theory and how their representation influences the shape of *g*(*r*) due to particle crowding and/or attraction. Comparisons between RPPs and KS-MD are done in Sec. III, where we begin by comparing interparticle forces illustrating how an increase in temperature indicates an increase in accuracy. In addition, the microfield distribution of forces, Einstein frequency, power spectrum, self-diffusion coefficient, and *g*(*r*) are compared, highlighting how an approximately accurate *g*(*r*) does not ensure similar accuracy in time correlation functions and transport coefficients. A description of how we accurately compute the self-diffusion coefficient and its uncertainty when finite-size errors are non-negligible is given in Sec. III C. This further emphasizes the need for RPPs, as we minimize finite-size errors in KS-MD simulations by making the necessary corrections as shown in Sec. III E. We conclude by comparing fully converged (in particle number and simulation time) self-diffusion coefficients to an analytic transport theory; benchmarking its accuracy and providing an effective interaction correction to extend the range of applicability.

## II. MODELS FOR THE INTERACTION POTENTIAL

The theoretical foundations of the models we will compare are described in this section; their connections are shown in Fig. 1. We compare three classes of interactions that are based on the ionic *N*-body energy, shown in the top box, pair interactions that are pre-computed and are analytic or tabulated, shown in the lower-left box, and optimal pair interactions extracted from the *N*-body results, shown in the lower-right box. By comparing these three approaches, we aim to answer several specific questions. First, given the nuclear charge *Z*, ionic number density *n _{i}*, and temperature

*T*, what ranges in ${Z,ni,T}$ space are the fast, pre-computed interactions valid and therefore allow for large-scale heterogeneous simulations? Second, how accurate is the “optimal” pair interaction, and what do its limitations reveal about the need for three-body interactions (and perhaps beyond)? Can these interactions be used to test and correct for finite-size errors? Third, can the optimal interactions guide the development of pre-computed interactions? To simplify the discussion we will consider single species matter with a range of

*Z*, each species at its normal solid ionic mass density

*ρ*, or in some cases half of that, and in thermodynamic equilibrium at temperature

_{i}*T*. While we do not consider mixtures in this work, the framework is general and can be straightforwardly applied to them.

Assuming the Born–Oppenheimer approximation holds, we define a potential energy surface for the ions as

Physically, the ionic potential energy surface is determined by the electronic charge distribution arising from ions at a particular set of coordinates; in general, (1) does not simplify into sums over pairwise terms. There are two major approaches to obtaining (1) in practice. The approach represented by the top box in Fig. 1 computes the electronic charge distribution for each ionic distribution. This is achieved computationally in Kohn–Sham approaches by reducing the electron many-body problem to a single-electron problem in which the Kohn–Sham electron moves in the external field of *N*-ionic centers. The dominant computational cost comes from solving an $No\xd7No$ set of eigenvalue equations, where *N _{o}* is the number of single-particle orbitals. Even though the electron many-body problem has been simplified to a one-body problem, matrix diagonalization incurs a cost of $O(No3)$, and at high temperatures the smearing of the Fermi–Dirac distribution requires an increasing number of orbitals leading to significant increases in computational cost. The complexity of the electron charge distribution also demands the use of an advanced “Jacob's ladder” of exchange-correlation (XC) functions to address the electron many body problem.

This approach yields an intrinsically ionic *N*-body potential energy surface; the electronic density is computed using a description appropriate to the choice of ${Z,ni,T}$. The second approach to calculating the potential energy surface is to use a cluster-type expansion, which takes the form

When this expansion can be truncated with only a few terms, interactions can be pre-computed, and fast neighbor algorithms allow for a very rapid evaluation of forces, typically many orders of magnitude faster than through use of (1). This allows, for example, for simulations with trillions of particles.^{16–18} However, the disadvantages are that the computational cost increases rapidly as more terms are included, and the accuracy of a specific truncation and choice of functional forms with that truncation are not usually known; part of our goal is to assess how accurate the potential energy surface in (1) can be represented by the first two terms of (2).

### A. *N*-body interaction potentials

The most accurate forces are obtained from the gradient of the total energy in (1), which requires the entire ionic configuration. Although machine learning approaches are enabling the ability to pre-learn that relationship,^{19–21} it remains more common to compute the forces for each ionic configuration during the simulation (“on-the-fly”). We obtain the electronic number density for each ionic configuration in the Kohn–Sham–Mermin formulation of the density

where *T* is the temperature of the system in energy units, the Fermi occupations are given by $fi(T)=(1+e\beta (Ei\u2212\mu ))\u22121$, and the Kohn–Sham–Mermin orbitals $\varphi i(r)$ satisfy

where

is a sum of the external (*N* ion–electron), Hartree, and exchange-correlation energies. Our KS-MD simulations were done using the Vienna *Ab initio* simulation package (VASP).^{22–25} The finite temperature electronic structure was treated with the Mermin free-energy functional, and we used the Perdew–Burke–Ernzerhof (PBE) functional for the exchange correlation energy.^{26} To improve computational efficiency, we eliminated the chemically inactive core electrons with the projector augmented-wave^{27} pseudopotential. Due to the anticipated high temperatures and small interionic separations, we used the smallest core “GW” pseudopotentials available in VASP. Sixty-four atoms (*N *=* *64) were used in these simulations, with an energy cutoff of 800 eV and at the Baldereschi mean-value k-point^{28} for all temperatures ranging from *T *=* *0.5 to 15 eV. A simulation time step of 0.1 fs was used, and the total simulation lengths for each case vary and are on the order of a few picoseconds. All KS-MD simulations were first equilibrated in the NVT ensemble and then carried out in the NVE ensemble where data were collected.

### B. Force matching

After the Kohn–Sham potential energy surface has been computed, we aim to construct a compact representation of (1) with (2). By assuming a parameterized functional form for (2), the force-matching procedure^{29–34} was used to generate the optimal RPP model based on the KS-MD force data. From each KS-MD simulation, a dataset of $K\u22613NM$ forces (3 force components, *N* atoms, and *M* atomic configurations) is obtained. Atomic simulation data at nearby time points are highly correlated; thus, a stride between atomic configurations was used to generate 100–200 independent configurations.

With each KS-MD dataset, we determine the optimal RPP for that system by minimizing the loss function

Here, *ζ* is a set of optimizable parameters, $Fk(\zeta )$ is the *k*th force for the parameterized model with parameters *ζ*, $Fk0$ is the *k*th force from KS-MD reference dataset, and *w _{k}* is a weight factor. The weight factor $wk=1/(Fk0+\epsilon )2$ ensures that both large and small forces contribute equally to (6). The parameter

*ε*should be varied for each temperature and element but in most cases here, $\epsilon \u22481$.

The choice of parameterization can either have a pre-computed functional form such as (8) or be determined completely from the data as is the case for a tabulated potential^{35} with spline interpolation—the choice in this work. For each system, we begin by sampling a TFY RPP at 15 locations in *r* and use that as the initial condition for the force-matching procedure. The TFY RPP is sampled such that $rmin<r<8\u2009\xc5$ where *r _{min}* is the minimum ionic separation in the KS-MD dataset. To ensure that the core repulsion and/or attractive oscillation regions are sampled sufficiently, 10 points are placed in the region where $rmin<r<4\u2009\xc5$, leaving the 5 remaining points to be placed where $r>4\u2009\xc5$. To test for convergence of the optimal force-matched RPP, two optimization methods were used (specifically simulated annealing and differential evolution). By choosing a tabulated potential form for the RPP, the explicit form of the model is entirely determined from the KS-MD force data and not limited to a fixed functional form.

While the force-matched RPP yields the best RPP to reproduce the KS-MD force data, it could be the case three-body and higher interactions are non-negligible. To check this, we selectively employ the spectral neighborhood analysis potential (SNAP) which constructs a potential energy surface from a set of four-body descriptors (bispectrum components), where each descriptor is independently weighted, and these weights are determined by regressing against KS-MD data of energies and forces. A descriptor captures the strength of density correlations between neighboring atoms and the central atom within a given cutoff distance; details can be found in Refs. 36 and 37. The parameterization of the SNAP uses descriptors of the local atomic environment capturing up to four-body interactions when represented in the form of (2), so lower errors associated with SNAP compared to an optimal RPP are entirely due to three- and four-body interactions. While higher bodied inter-atomic potentials exist in the literature,^{38} it can be expected there are diminishing accuracy returns with higher interaction moments; thus, SNAP offers a leading order check on the RPP compared here.

SNAP potentials utilizing 56 bispectrum component descriptors were trained on 10% of the KS-MD dataset and additionally tested against an additional 10% to ensure regression errors were properly minimized and avoided over-fitting of the KS-MD data.

### C. Radial pair potentials

As the computational cost of using on-the-fly *N*-body interactions is often prohibitive, the least expensive approach utilizes pre-computed RPPs ignoring most of the terms in (2). Many functional forms for the RPP have been proposed for application to warm dense matter often using the second-order perturbation-theory interaction energy

which is the standard Fourier-space result^{39} written in terms of the mean ionization state $\u27e8Z\u27e9$, the bare Coulomb potential $uC(k)=4\pi e2/k2$, the electron–ion pseudopotential $uei(k)$, and the susceptibility $\chi (k)$.

In practice, pair interactions are constructed using nearly the same steps as for the *N*-body interactions, with the primary difference being that each ion is replaced with a single “average atom” (AA), which is an all-electron, non-linear, finite-temperature density functional theory calculation;^{40} such calculations can also be relativistic.^{41,42} From the AA, a pseudopotential $uei(k)$ and an accurate free/valence electron response function $\chi (k)$ are constructed and (7) is formed. This approach has three strengths: (1) typical AA models are not limited to low temperatures, (2) the interaction (7) can be pre-computed for use in MD, and (3) pair interactions with a fast nearest neighbor algorithm are very computationally efficient. As we alluded to above, the accuracy loss attendant to these strengths is what we wish to determine in this work. The AA itself is aware of the ionic number density *n _{i}*, which sets the ion-sphere radius $ai=(3/4\pi ni)1/3$, and includes the fact that there is only one ion in the ion sphere, which implies a

*g*(

*r*); this indirect inclusion of higher-order terms in (2) is true for all AA-based interactions.

Among the simplest variants of (7), one approximates the pseudopotential as $uei(k)\u2248\u22124\pi \u27e8Z\u27e9e2/k2$, where the mean ionization state $\u27e8Z\u27e9$ results from a AA calculation,^{40} and $\chi (k)$ in its long-wavelength (Thomas–Fermi) limit $\chi TF(k)$; this is known as the “Yukawa” interaction.^{10,43} Here, we employ a Yukawa interaction with inputs from a Thomas–Fermi AA,^{1} which we will refer to as “TFY.” This procedure yields an analytic potential in real space of the form

where the electron screening is approximated by the Thomas–Fermi screening length

where $F\u22121/2$ is the Fermi–Dirac integral of order $\u22121/2,\u2009\beta =1/T$, and *μ _{e}* is the electron chemical potential. Padé fits of Fermi–Dirac integrals and their inverses are carried out in Refs. 44 and 45. An approximation to these fits

^{46}yields

where the Fermi energy $EF=\u210f2(3\pi 2ne)2/3/2me$. Note that the TFY interaction is monotonically decreasing (purely repulsive). Computationally, the TFY model is highly desirable because of its radial, pair, analytic form with an exponentially damped short range. Its weaknesses are the relatively approximate treatments of $uei(k)$ and $\chi (k)$. The TFY model can be extended by including the gradient corrections to $\chi TF(k)$, but otherwise retaining the other approximations. This improvement yields the Stanton–Murillo potential;^{10} the gradient correction to $\chi TF(k)$ introduces oscillations in the potential in some plasma regimes that are absent in the monotonic TFY model. Moreover, gradient corrections add improvements to the cusp at the origin and the large-*r* asymptotic behavior. Here, however, we will only employ the simpler TFY model.

A great deal of accuracy can be gained by abandoning analytic inputs to (7). In this case, self-consistent numerical calculations of each of the terms can be carried out, still allowing for pre-computed interactions; there is essentially no computational overhead for tabulated interactions.^{35} Here, we employ a NPA model that yields both the mean ionization state and its pseudopotential using a Kohn–Sham–Mermin approach, as described above, but with a finite-temperature exchange-correlation potential; the susceptibility is determined by the Lindhard function with local field corrections.^{13} Note that the electron–ion pseudopotential $uei(k)$ introduces additional oscillations on length scales different from $\chi (k)$ although the Friedel oscillations in $\chi (k)$ contribute much more to the pair interaction. Note that the name “NPA” has been used by many authors to several different average-atom models, and many of them involve approximations that limit those models to higher temperatures, e.g., $T>EF$; however, here we use the one-center density functional theory model developed by Dharma-wardana and Perrot as this model has been tested at high temperatures as well as at very low temperatures, and found to agree closely with more detailed *N*-center density functional theory simulations and path-integral quantum calculations where available.

It is worth comparing predictions based on (7) with other forms suggested previously. A popular RPP for warm dense matter studies is the short-range repulsion interaction, which adds a long-range, power-law correction to the TFY model of the form $A/r4$;^{9,47–52} for *A *>* *0, this is also a monotonic interaction, with the goal of increasing the strength of the TFY model, which underestimates the peak height of *g*(*r*). In Fig. 2, we examine this ansatz by computing a NPA interaction for Al at solid density and *T *=* *1 eV. To find the “best” power law, we multiply the NPA interaction by various powers *r ^{a}* to find regions where the interaction is flat; a flat region with

*a*=

*4 would recover the short-range repulsion interaction. It is clear that the $A/r4$ is only valid over a very small range of*

*r*values; importantly, the NPA interaction shows that the exponent

*a*increases as

*r*becomes large, which is a true short-ranged interaction—the empirical correction the short-range repulsion model adds greatly overestimates the strength of the interaction at large interparticle separations.

^{12}Worse, the short-range repulsion model potentially gets an accurate answer for the wrong reason, as we explore in Fig. 3.

Because the form (7) generally has oscillations, the enhanced peak height of *g*(*r*) from the NPA model over the TFY model occurs for two, independent reasons. Attractive regions of the interaction, as shown in the top panel of Fig. 3, can produce very strong peaks in *g*(*r*). Conversely, stronger overall repulsion at intermediate *r* can lead to a similar *g*(*r*) behavior, as shown in the bottom panel of Fig. 3, but with rapid decay of the interaction at larger *r*. The functional form (7) naturally contains both the “crowding” and “attraction” behaviors as special cases. Figure 4 shows a comparison of the RPPs for C, Al, V, and Au at *T *=* *0.5 and 5 eV. The TFY model is purely monotonic, whereas the force-matched and NPA RPPs have attractive and repulsive regions in their oscillations. Below, we will explore the consequences of these features of the interaction on ionic transport.

Once the RPPs have been constructed, MD simulations were carried out using in the large-scale atomic/molecular massively parallel simulator (LAMMPS).^{53} For the tabulated RPPs (force-matched and NPA) a linear interpolation was needed to determine the force value between tabulation points. To make a direct comparison between the RPP-MD and KS-MD results, all simulations were carried out in a three dimensional periodic box with 64 atoms and a time step of 0.1 fs. The length of each simulation is identical to the corresponding simulation performed with KS-MD. Keeping these conditions identical avoids the unintentional reduction in statistical errors between KS-MD and RPP-MD. All simulations were first equilibrated in the NVT ensemble so that the average temperature for each simulation during the data collection phase is within 1% of the reported temperature in Table I. The data collection phase was carried out in the NVE ensemble. In Sec. III E, a finite-size effect study was done for the cases of C at 2.267 g/cm^{3} and V at 6.11 g/cm^{3} where the total simulation length was increased by 10 times and the number of atoms *N* increases from 64 to 256, 3375, and 8000.

Element . | ρ (g/cm_{i}^{3})
. | T (eV)
. | $DKS\u2212MD\u2009(10\u22123\xc52/fs)$ . | $DFM\u2009(10\u22123\xc52/fs)$ . | $DTFY\u2009(10\u22123\xc52/fs)$ . | $DNPA\u2009(10\u22123\xc52/fs)$ . |
---|---|---|---|---|---|---|

Li | 0.513 | 0.054 | 1.4 ± 0.13 | 1.27 ± 0.054 | 5.6 ± 0.39 | 1.26 ± 0.077 |

C | 2.267 | 0.47 | 2.4 ± 0.12 | 9.3 ± 0.20 | 11.0 ± 0.60 | 1.69 ± 0.060 |

1.0 | 18.6 ± 0.7 | 27.2 ± 0.77 | 25 ± 1.64 | 11.0 ± 0.45 | ||

2.0 | 46 ± 1.68 | 49 ± 1.0 | 42 ± 3.70 | 43 ± 3.86 | ||

4.9 | 85 ± 5 | 94 ± 5.82 | 106 ± 5.37 | 92 ± 3.45 | ||

10.0 | 215 ± 4.49 | 151 ± 4.63 | ||||

15 | 266 ± 3.46 | 198 ± 5.10 | ||||

20 | 349 ± 7.60 | 249 ± 2.35 | ||||

28 | 423 ± 13.28 | 324 ± 12.17 | ||||

Al | 2.7 | 0.1 | 1.6 ± 0.14 | 0.35 ± 0.021 7 | ||

0.50 | 3.8 ± 0.16 | 4.1 ± 0.13 | 9.17 ± 0.099 | 3.9 ± 0.11 | ||

1.1 | 9.8 ± 0.30 | 9.4 ± 0.11 | 17.5 ± 0.70 | 8.5 ± 0.44 | ||

2.0 | 18.7 ± 0.50 | 18.8 ± 0.68 | 34.8 ± 0.52 | 18.8 ± 0.40 | ||

4.9 | 48 ± 3.56 | 49 ± 3.17 | 72 ± 3.03 | 54 ± 2.7 | ||

9.2 | 83 ± 1.63 | 84 ± 5.67 | 122 ± 5.83 | 94 ± 8.80 | ||

10.0 | 131 ± 6.77 | 105 ± 13.53 | ||||

15.4 | 134 ± 3.68 | 129 ± 3.37 | 169 ± 5.26 | 142 ± 6.17 | ||

20.0 | 197 ± 5.21 | 151 ± 2.55 | ||||

30.0 | 252 ± 4.74 | 203 ± 6.44 | ||||

Ar | 1.395 | 0.48 | 10.7 ± 0.43 | 12 ± 1.03 | 19 ± 1.10 | |

1.0 | 20.1 ± 0.89 | 26 ± 3.0 | 39 ± 2.22 | |||

2.0 | 48 ± 1.75 | 45 ± 2.84 | 85 ± 8.75 | |||

5 | 143 ± 6.55 | 171 ± 4.09 | ||||

10.0 | 210 ± 14.53 | 179 ± 6.21 | ||||

15.0 | 235 ± 13.34 | 193 ± 11.95 | ||||

20.0 | 255 ± 6.73 | 209 ± 8.91 | ||||

30.0 | 268 ± 2.73 | 228 ± 8.26 | ||||

V | 6.11 | 0.49 | 2.25 ± 0.050 | 2.86 ± 0.079 | 3.9 ± 0.18 | 0.91 ± 0.027 |

1.0 | 5.5 ± 0.21 | 6.5 ± 0.16 | 7.9 ± 0.36 | 6.6 ± 0.15 | ||

2.1 | 11.6 ± 0.78 | 12.5 ± 0.68 | 17.8 ± 0.74 | 14.8 ± 0.50 | ||

4.8 | 24.2 ± 0.63 | 24.7 ± 0.88 | 41 ± 2.76 | 27.7 ± 0.90 | ||

9.5 | 46 ± 3.41 | 42 ± 2.65 | 68 ± 2.10 | 47.6 ± 0.93 | ||

14.6 | 53 ± 1.81 | 57 ± 3.25 | 84 ± 4.83 | 63 ± 1.19 | ||

20.0 | 103 ± 6.10 | 82.7 ± 0.78 | ||||

30.0 | 134 ± 8.57 | 96 ± 1.86 | ||||

3.055 | 0.5 | 9.0 ± 0.81 | 11.3 ± 0.29 | 8.7 ± 0.23 | ||

0.97 | 14.7 ± 0.47 | 15.4 ± 0.43 | 19 ± 1.39 | |||

2.0 | 23 ± 1.13 | 24 ± 1.84 | 31 ± 1.26 | 27 ± 1.84 | ||

4.9 | 47 ± 4.38 | 44 ± 2.52 | 66 ± 7.30 | 48 ± 2.02 | ||

Fe | 7.874 | 0.51 | 2.13 ± 0.047 | 2.34 ± 0.042 | 2.84 ± 0.030 | |

1.1 | 5.27 ± 0.098 | 5.5 ± 0.16 | 5.9 ± 0.39 | |||

2.1 | 10.4 ± 0.72 | 10.4 ± 0.73 | 14.8 ± 0.46 | 9.2 ± 0.47 | ||

5.0 | 20.4 ± 0.61 | 22.0 ± 0.97 | 32 ± 1.41 | 27.1 ± 0.19 | ||

10.4 | 35 ± 1.14 | 38 ± 1.40 | 54 ± 1.25 | 49 ± 2.90 | ||

15.0 | 83 ± 5.18 | 60 ± 2.60 | ||||

20.0 | 97 ± 2.93 | 70 ± 1.04 | ||||

30.0 | 103 ± 4.69 | 83.0 ± 0.94 | ||||

3.937 | 0.51 | 6.0 ± 0.39 | 8.5 ± 0.94 | 6.2 ± 0.21 | ||

1.1 | 15.8 ± 0.70 | 15.6 ± 0.67 | 14.4 ± 0.33 | |||

2.1 | 20 ± 1.18 | 22 ± 2.07 | 27 ± 1.28 | |||

Au | 19.30 | 0.52 | 0.92 ± 0.028 | 0.71 ± 0.084 | 1.67 ± 0.12 | 0.51 ± 0.042 |

1.1 | 2.0 ± 0.11 | 1.92 ± 0.088 | 3.9 ± 0.42 | 1.66 ± 0.069 | ||

1.9 | 4.0 ± 0.14 | 3.4 ± 0.15 | 6.6 ± 0.16 | 3.52 ± 0.05 | ||

5.0 | 7.8 ± 0.40 | 8.2 ± 0.21 | 15.3 ± 0.63 | 10.7 ± 0.50 | ||

9.7 | 14.4 ± 0.64 | 15 ± 1.19 | 25 ± 1.94 | 16.6 ± 0.86 | ||

15.0 | 19.82 ± 0.80 | 22 ± 2.56 | 30 ± 1.95 | 25 ± 2.79 | ||

20.0 | 39 ± 3.49 | 28.0 ± 0.97 | ||||

30.0 | 56 ± 2.23 | 33 ± 1.26 | ||||

Element . | ρ (g/cm_{i}^{3})
. | T (eV)
. | $DKS\u2212MD\u2009(10\u22123\xc52/fs)$ . | $DFM\u2009(10\u22123\xc52/fs)$ . | $DTFY\u2009(10\u22123\xc52/fs)$ . | $DNPA\u2009(10\u22123\xc52/fs)$ . |
---|---|---|---|---|---|---|

Li | 0.513 | 0.054 | 1.4 ± 0.13 | 1.27 ± 0.054 | 5.6 ± 0.39 | 1.26 ± 0.077 |

C | 2.267 | 0.47 | 2.4 ± 0.12 | 9.3 ± 0.20 | 11.0 ± 0.60 | 1.69 ± 0.060 |

1.0 | 18.6 ± 0.7 | 27.2 ± 0.77 | 25 ± 1.64 | 11.0 ± 0.45 | ||

2.0 | 46 ± 1.68 | 49 ± 1.0 | 42 ± 3.70 | 43 ± 3.86 | ||

4.9 | 85 ± 5 | 94 ± 5.82 | 106 ± 5.37 | 92 ± 3.45 | ||

10.0 | 215 ± 4.49 | 151 ± 4.63 | ||||

15 | 266 ± 3.46 | 198 ± 5.10 | ||||

20 | 349 ± 7.60 | 249 ± 2.35 | ||||

28 | 423 ± 13.28 | 324 ± 12.17 | ||||

Al | 2.7 | 0.1 | 1.6 ± 0.14 | 0.35 ± 0.021 7 | ||

0.50 | 3.8 ± 0.16 | 4.1 ± 0.13 | 9.17 ± 0.099 | 3.9 ± 0.11 | ||

1.1 | 9.8 ± 0.30 | 9.4 ± 0.11 | 17.5 ± 0.70 | 8.5 ± 0.44 | ||

2.0 | 18.7 ± 0.50 | 18.8 ± 0.68 | 34.8 ± 0.52 | 18.8 ± 0.40 | ||

4.9 | 48 ± 3.56 | 49 ± 3.17 | 72 ± 3.03 | 54 ± 2.7 | ||

9.2 | 83 ± 1.63 | 84 ± 5.67 | 122 ± 5.83 | 94 ± 8.80 | ||

10.0 | 131 ± 6.77 | 105 ± 13.53 | ||||

15.4 | 134 ± 3.68 | 129 ± 3.37 | 169 ± 5.26 | 142 ± 6.17 | ||

20.0 | 197 ± 5.21 | 151 ± 2.55 | ||||

30.0 | 252 ± 4.74 | 203 ± 6.44 | ||||

Ar | 1.395 | 0.48 | 10.7 ± 0.43 | 12 ± 1.03 | 19 ± 1.10 | |

1.0 | 20.1 ± 0.89 | 26 ± 3.0 | 39 ± 2.22 | |||

2.0 | 48 ± 1.75 | 45 ± 2.84 | 85 ± 8.75 | |||

5 | 143 ± 6.55 | 171 ± 4.09 | ||||

10.0 | 210 ± 14.53 | 179 ± 6.21 | ||||

15.0 | 235 ± 13.34 | 193 ± 11.95 | ||||

20.0 | 255 ± 6.73 | 209 ± 8.91 | ||||

30.0 | 268 ± 2.73 | 228 ± 8.26 | ||||

V | 6.11 | 0.49 | 2.25 ± 0.050 | 2.86 ± 0.079 | 3.9 ± 0.18 | 0.91 ± 0.027 |

1.0 | 5.5 ± 0.21 | 6.5 ± 0.16 | 7.9 ± 0.36 | 6.6 ± 0.15 | ||

2.1 | 11.6 ± 0.78 | 12.5 ± 0.68 | 17.8 ± 0.74 | 14.8 ± 0.50 | ||

4.8 | 24.2 ± 0.63 | 24.7 ± 0.88 | 41 ± 2.76 | 27.7 ± 0.90 | ||

9.5 | 46 ± 3.41 | 42 ± 2.65 | 68 ± 2.10 | 47.6 ± 0.93 | ||

14.6 | 53 ± 1.81 | 57 ± 3.25 | 84 ± 4.83 | 63 ± 1.19 | ||

20.0 | 103 ± 6.10 | 82.7 ± 0.78 | ||||

30.0 | 134 ± 8.57 | 96 ± 1.86 | ||||

3.055 | 0.5 | 9.0 ± 0.81 | 11.3 ± 0.29 | 8.7 ± 0.23 | ||

0.97 | 14.7 ± 0.47 | 15.4 ± 0.43 | 19 ± 1.39 | |||

2.0 | 23 ± 1.13 | 24 ± 1.84 | 31 ± 1.26 | 27 ± 1.84 | ||

4.9 | 47 ± 4.38 | 44 ± 2.52 | 66 ± 7.30 | 48 ± 2.02 | ||

Fe | 7.874 | 0.51 | 2.13 ± 0.047 | 2.34 ± 0.042 | 2.84 ± 0.030 | |

1.1 | 5.27 ± 0.098 | 5.5 ± 0.16 | 5.9 ± 0.39 | |||

2.1 | 10.4 ± 0.72 | 10.4 ± 0.73 | 14.8 ± 0.46 | 9.2 ± 0.47 | ||

5.0 | 20.4 ± 0.61 | 22.0 ± 0.97 | 32 ± 1.41 | 27.1 ± 0.19 | ||

10.4 | 35 ± 1.14 | 38 ± 1.40 | 54 ± 1.25 | 49 ± 2.90 | ||

15.0 | 83 ± 5.18 | 60 ± 2.60 | ||||

20.0 | 97 ± 2.93 | 70 ± 1.04 | ||||

30.0 | 103 ± 4.69 | 83.0 ± 0.94 | ||||

3.937 | 0.51 | 6.0 ± 0.39 | 8.5 ± 0.94 | 6.2 ± 0.21 | ||

1.1 | 15.8 ± 0.70 | 15.6 ± 0.67 | 14.4 ± 0.33 | |||

2.1 | 20 ± 1.18 | 22 ± 2.07 | 27 ± 1.28 | |||

Au | 19.30 | 0.52 | 0.92 ± 0.028 | 0.71 ± 0.084 | 1.67 ± 0.12 | 0.51 ± 0.042 |

1.1 | 2.0 ± 0.11 | 1.92 ± 0.088 | 3.9 ± 0.42 | 1.66 ± 0.069 | ||

1.9 | 4.0 ± 0.14 | 3.4 ± 0.15 | 6.6 ± 0.16 | 3.52 ± 0.05 | ||

5.0 | 7.8 ± 0.40 | 8.2 ± 0.21 | 15.3 ± 0.63 | 10.7 ± 0.50 | ||

9.7 | 14.4 ± 0.64 | 15 ± 1.19 | 25 ± 1.94 | 16.6 ± 0.86 | ||

15.0 | 19.82 ± 0.80 | 22 ± 2.56 | 30 ± 1.95 | 25 ± 2.79 | ||

20.0 | 39 ± 3.49 | 28.0 ± 0.97 | ||||

30.0 | 56 ± 2.23 | 33 ± 1.26 | ||||

## III. RESULTS

### A. Force error analysis

One metric for establishing the accuracy of approximations to the Kohn–Sham potential energy surface is to compute relative force errors between Kohn–Sham force data and a parameterized model (RPP or many-body potential) for *M* particle coordinate configurations. For this, we compute the mean-absolute force error

where $F\alpha ,i,m(PAR)$ and $F\alpha ,i,m(KS\u2212MD)$ are the *α*th force components (*x*, *y*, or *z*) on the *i*th atom in particle coordinate configuration number *m* for the parameterized model and the KS-MD force data, respectively.

Note that a direct comparison of the mean absolute error between different elements, temperatures, and densities cannot be done as the distributions of forces associated with systems of different elements at different thermodynamic conditions are in general quite different. This can be observed in Fig. 5 where a microfield distribution of the force magnitudes is shown. In all cases but C at 2.267 g/cm^{3} and *T *=* *5 eV, the TFY model peaks at a smaller field value than KS-MD. In contrast, for C, V, and Au at *T *=* *0.5 eV, the NPA RPP peaks at a higher field value than KS-MD. These trends can be connected back to (7) where the choice of $\u27e8Z\u27e9,\u2009uei(k)$, and $\chi (k)$ all contribute to the construction of a RPP model and hence the force magnitudes. More work needs to be done to determine how each term influences the RPP model, the predicted forces, and observables.

As the microfield force distributions vary for different elements and temperatures, the mean absolute error will also vary. To this end, we seek a scale factor for (11) to normalize the results across the different elements, temperatures, and densities studied here. Such a scale factor is the “mean absolute force” defined as

This metric has the following desirable property: if the mean absolute error changes with density or temperature in the same way as the underlying force distribution, the relative force error will maintain roughly the same value. Therefore, as we change the thermodynamic conditions for a given element, (13) provides a temperature independent metric as measured with respect to a KS-MD force data “baseline.” Intuitively, when (13) evaluates to 1, the mean absolute error is the same order of magnitude as the mean absolute force and when $(13)$ is zero, the parameterized model is exactly reproducing the per-component KS-MD force data.

Figure 6 displays (13) as a function of temperature for C, Al, V, and Au where general trends can be observed. One trend is that for most RPPs, the relative force error decreases toward higher temperatures, which confirms an intuition long held for the validity of the NPA and TFY models. However, for all systems pictured except C, force matching drastically reduces the relative force error compared to the NPA and TFY results. Specifically, the force-matched RPPs routinely achieve a relative force error of roughly 0.05 above *T *=* *5 eV. Except for the case of the NPA RPP for Al, the NPA and TFY RPPs maintain an error of around 0.2 across the entire the temperature range.

The second major observation from Fig. 6 is that while force-matched RPPs drastically lower the observed relative force errors across temperatures compared against other RPPs, we immediately see where a RPP approximation is likely invalid. For example, the relative force error for C using the force-matched RPP is uncharacteristically high (roughly 0.6) until *T *=* *5 eV. A similar situation appears for the case of V at *T *=* *0.5 eV where the relative force error for the force-matched RPP is roughly 0.25. We can demonstrate explicitly that these discrepancies come from the neglect of three-body and higher interactions by showing relative force errors using a SNAP model. For C, the relative force error drops from roughly 0.6 using a force-matched RPP to 0.2 using a SNAP model at *T *=* *0.5 eV. Likewise for V, the relative force error drops from roughly 0.25 using a force-matched RPP to 0.07 using a SNAP model at the same temperature.

Ultimately, it is not the component-wise force or the interaction potential we care about generating, but rather observables such as *g*(*r*) and the self-diffusion coefficient. To address this connection, we examine correlations between the force error and the self-diffusion coefficient error, as shown in Fig. 7. While there is a general trend with increasing errors in both quantities (shown with a linear fit), there are also some clear outliers. For the case of C at 2.267 g/cm^{3} and *T *=* *0.5 eV, we find that the NPA and TFY RPPs produce a self-diffusion coefficient that differs from the KS-MD result by many factors. However, C under these conditions exists in several charge states with transient bonding; only the NPA accounts for this. (More details are found in the Appendix.) This case is marked with arrows in Fig. 7. Conversely, for V at *T *=* *1 eV, the relative self-diffusion percent error is low, yet the relative forc e error is high. The imperfect mapping of relative self-diffusion error vs relative force error suggests that physics beyond a RPP is needed, possibly at least a three-body angular dependence, but further work is needed.

### B. Radial distribution function and the Einstein frequency

The radial distribution function^{30} is a measure of spatial correlations normalized by the ideal gas. It has been shown that, in general, there always exists a RPP that can reproduce *g*(*r*) from a *N*-body simulation,^{54} and the force-matching procedure provides an avenue for obtaining this RPP. Figure 8 compares *g*(*r*) computed from MD simulations for all RPP models for C, Al, V, and Au. Each row corresponds to a different temperature, and clear trends can be observed, such as the improvement in agreement between models as the temperature increases. We note that the force-matched RPP always obtains the correct *g*(*r*), and the NPA model generally predicts the location of the first peak but sometimes over-predicts the magnitude or misses the location of the first peak altogether as observed in the case of V at 6.11 g/cm^{3} for *T *=* *0.5 eV. The TFY model always underestimates the magnitude of the first peak height, and the location is usually shifted.

Insight into the connection between the *g*(*r*) peak height and the self-diffusion coefficient can be obtained from the normalized velocity autocorrelation function^{55}

where $v(t)$ is the velocity of a particle at time *t* and $\u27e8\xb7\u27e9$ is an ensemble average over particles *and* time. A short time expansion of (14) yields

where $\Omega 0$ is the Einstein frequency

where *m _{i}* is the ion mass in grams. The Einstein frequency gives insight into the relationship between

*u*(

*r*) and

*g*(

*r*), highlighting how different regions are weighted more or less depending on the curvature of

*u*(

*r*). In Fig. 9, the integrand of (16) is shown. For the TFY model, the integrand is always smaller than those predicted by force-matched and NPA RPPs. The area under each curve in Fig. 9 can be directly connected to the self-diffusion coefficient through the Green–Kubo relation (in three dimensions)

substituting (15) into (17). Doing so shows that the TFY model will always predict a larger self-diffusion coefficient than the force-matched or NPA model as the area under these curves is larger. This is confirmed later when the self-diffusion coefficients are explicitly calculated as discussed in Sec. III C.

### C. Self-diffusion

Another approach to compute the self-diffusion coefficient is via the slope of the mean-squared displacement from the Einstein relation

Both (17) and (18) can be used to compute the self-diffusion coefficient and have been shown to be equivalent.^{56} In this work, the self-diffusion coefficient has been calculated from a linear fit to the mean-squared displacement, $\u27e8|r(t)\u2212r(0)|2\u27e9$.

Due to finite-size effects, two problems arise when computing the slope and uncertainty of the linear fit. First, we must ensure that the linear fit is carried out in the late-time linear regime of the mean-squared displacement. Second, we dismiss statistically unconverged late-time behavior of the mean-squared displacement where the ensemble average contains sparse amounts of data. To remedy both of these concerns, we uniformly randomly sub-sample the mean-squared displacement 100 times with 10 points along each sub-sample. Next, a linear fit is determined for each sub-sample, and the standard deviation of the sub-sample slopes is computed. Once the standard deviation is known, a cutoff time is calculated by determining the point in time that the standard deviation of the sub-sample fits is less than half of the standard deviation computed from *sub-sample* fits to the *entire* mean-squared displacement. The simulation data for the mean-squared displacement after the cutoff time is discarded, and the fitting procedure described above is repeated. The average, and standard deviation of the fits to the reduced dataset yield self-diffusion coefficient and the uncertainty, respectively, and are reported in Table I and displayed in Fig. 10.

Given the values for the self-diffusion coefficient reported in Table I, we can answer the following question: at what temperature are computationally inexpensive models adequate? To do this, we compute the relative self-diffusion coefficients $DNPA/DKS\u2212MD$ and $DTFY/DNPA$. For example, the top panel in Fig. 11 suggests that NPA models may be accurate from *T *=* *1 eV and above if the target error tolerance is 50% of the self-diffusion coefficient computed from KS-MD. Similarly in the bottom figure, the TFY model is generally accurate to within 50% of the NPA model from *T *=* *5 eV and beyond.

Two important observations can be made from the trends in Fig. 11. The top panel illustrates temperatures at which an *N*-body potential is needed and when NPA is adequate. The bottom panel shows a comparison with TFY, which has the simplest $uei(k)$ and $\chi (k)$, and we see temperatures at which TFY becomes comparable to NPA, suggesting when we can exploit simpler approximations for those inputs.

### D. Power spectrum

The self-diffusion coefficient is useful for comparing and quantifying the accuracy of RPP models and transport theories, but in order to assess how accurately the particle dynamics are reproduced, we look at the power spectrum of the velocity autocorrelation function *Z*(*t*)

In Fig. 12, we compare $Z\u0303(\nu )$ calculated using TFY, force-matched, and NPA RPPs against results obtained from KS-MD. We find that with the exception of low temperature C and V, force-matched RPPs agree with the KS-MD results across the entire frequency range. This, combined with the low relative force errors and accurate reproduction of static properties discussed previously, indicates that the force-matched RPPs accurately approximate the Kohn–Sham potential energy surface. For higher temperatures, the NPA RPP is very similar to the force-matched RPP for low and high frequencies for all elements. For *T *=* *0.5 eV, the dynamics predicted from the NPA model are noticeably less similar to those from KS-MD where NPA underestimates the prevalence of low-frequency modes in Au and both low and high-frequency modes in V. Interestingly, the NPA RPP captures the single-particle dynamics of low temperature C very well, but Figs. 4 and 8 indicate that this agreement comes at the expense of sacrificing the accuracy of static properties. Finally, the TFY RPP exhibits roughly the same trends across all elements and temperatures—overestimation of the low frequency modes and underestimation of the high-frequency modes except for the case of C at 2.267 g/cm^{3} and *T *=* *5 eV where excellent agreement with KS-MD is observed.

### E. Finite-size corrections

Generally, thousands or even millions of atoms are needed to approximate the thermodynamic limit.^{1,16} While the KS-MD framework provides an accurate description of the electronic structure and the *N*-body potential is determined on-the-fly, corrections for finite-size effects must be considered. When the shear viscosity *η* of the system is known, finite-size corrections can be determined from Ref. 57

where $D\u221e$ is the self-diffusion coefficient in the thermodynamic limit, *D _{N}* is the self-diffusion coefficient computed from a system of finite number of particles

*N*, and $\xi =2.837\u2009297$ for cubic simulation boxes of length L with periodic boundary conditions. When

*η*is unknown, multiple simulations of increasing particle number are carried out, and a linear fit is used to determine $D\u221e$. Results from this procedure are shown in Fig. 13 where $D\u221e$ is determined via linear extrapolation to $1/L=0$.

By finding the percent difference in $D\u221e$ and *D _{N}*, we approximate the errors from finite-size effects in the KS-MD self-diffusion coefficient at these conditions. The approximate error in KS-MD for the case shown in Fig. 13, is $\u223c20$%. While the error will vary with {

*Z*,

*n*,

*T*}, the impact of finite-size effects is significant. From this study, the most promising approach is to fully converge the NPA MD results, using force-matched RPPs when necessary (for low temperatures $T\u2009\u2272\u20091$ eV).

Finite-size corrections allow for a direct comparison to analytic transport theories, namely, the Stanton–Murillo model.^{58} The Stanton–Murillo model provides a closed form solution for ionic self-diffusion by using an effective interaction potential in a Boltzmann kinetic theory framework. The major benefit of this model is that the computation of ionic transport is nearly instantaneous. However, its applicability in the cold dense matter and warm dense matter regimes is unknown.

The results in Table II show that the effective interaction approach of the Stanton–Murillo model captures much of the many-body physics included in the TFY RPP results. The main weakness of the model, and also TFY, is therefore the functional form of the interaction they employ as the differences with the force-matched and NPA columns reveal. Because self-diffusion is a relatively simple transport coefficient,^{58} more work is needed to quantify these trends for other transport properties.

Element . | T (eV)
. | $DFM\u2009(10\u22123\xc52/fs)$ . | $DNPA\u2009(10\u22123\xc52/fs)$ . | $DTFY\u2009(10\u22123\xc52/fs)$ . | $DSM\u2009(10\u22123\xc52/fs)$ . | $DCSM\u2009(10\u22123\xc52/fs)$ . |
---|---|---|---|---|---|---|

C | 0.47 | 10.55 | 2.14 | 12.66 | 13.08 | 2.14 |

1.0 | 32.44 | 14.21 | 25.57 | 26.11 | 13.87 | |

2.0 | 56.70 | 43.12 | 51.14 | 50.53 | 39.23 | |

4.9 | 99.51 | 109.55 | 117.88 | 118.34 | 106.76 | |

10.0 | 169.91 | 210.76 | 217.34 | 206.71 | ||

15.0 | 219.99 | 296.21 | 293.54 | 284.10 | ||

20.0 | 256.33 | 342.15 | 356.41 | 348.04 | ||

28.0 | 327.32 | 470.10 | 439.44 | 432.07 | ||

V | 0.49 | 4.14 | 1.01 | 5.42 | 6.76 | 4.39 |

1.0 | 8.54 | 8.53 | 11.53 | 12.26 | 7.96 | |

2.1 | 15.67 | 18.87 | 22.56 | 23.14 | 15.03 | |

4.8 | 28.72 | 31.34 | 42.42 | 46.25 | 30.18 | |

9.5 | 49.49 | 54.90 | 73.76 | 77.10 | 51.16 | |

14.6 | 66.98 | 74.44 | 99.82 | 99.78 | 67.97 | |

20.0 | 87.90 | 118.24 | 117.8 | 83.15 | ||

30.0 | 105.60 | 143.63 | 141.66 | 106.94 | ||

50.0 | 131.84 | 175.35 | 171.07 | 143.02 | ||

75.0 | 178.34 | 202.93 | 194.31 | 172.91 | ||

100.0 | 209.50 | 207.20 | 211.86 | 194.36 |

Element . | T (eV)
. | $DFM\u2009(10\u22123\xc52/fs)$ . | $DNPA\u2009(10\u22123\xc52/fs)$ . | $DTFY\u2009(10\u22123\xc52/fs)$ . | $DSM\u2009(10\u22123\xc52/fs)$ . | $DCSM\u2009(10\u22123\xc52/fs)$ . |
---|---|---|---|---|---|---|

C | 0.47 | 10.55 | 2.14 | 12.66 | 13.08 | 2.14 |

1.0 | 32.44 | 14.21 | 25.57 | 26.11 | 13.87 | |

2.0 | 56.70 | 43.12 | 51.14 | 50.53 | 39.23 | |

4.9 | 99.51 | 109.55 | 117.88 | 118.34 | 106.76 | |

10.0 | 169.91 | 210.76 | 217.34 | 206.71 | ||

15.0 | 219.99 | 296.21 | 293.54 | 284.10 | ||

20.0 | 256.33 | 342.15 | 356.41 | 348.04 | ||

28.0 | 327.32 | 470.10 | 439.44 | 432.07 | ||

V | 0.49 | 4.14 | 1.01 | 5.42 | 6.76 | 4.39 |

1.0 | 8.54 | 8.53 | 11.53 | 12.26 | 7.96 | |

2.1 | 15.67 | 18.87 | 22.56 | 23.14 | 15.03 | |

4.8 | 28.72 | 31.34 | 42.42 | 46.25 | 30.18 | |

9.5 | 49.49 | 54.90 | 73.76 | 77.10 | 51.16 | |

14.6 | 66.98 | 74.44 | 99.82 | 99.78 | 67.97 | |

20.0 | 87.90 | 118.24 | 117.8 | 83.15 | ||

30.0 | 105.60 | 143.63 | 141.66 | 106.94 | ||

50.0 | 131.84 | 175.35 | 171.07 | 143.02 | ||

75.0 | 178.34 | 202.93 | 194.31 | 172.91 | ||

100.0 | 209.50 | 207.20 | 211.86 | 194.36 |

With the converged self-diffusion data, we generate an effective interaction correction to the Stanton–Murillo model. The effective interaction corrected Stanton–Murillo model is

where $\alpha (Z,T)$ is determined by fitting the ratio of the self-diffusion coefficient from the best performing RPP model and the self-diffusion coefficient computed from the Stanton–Murillo model *D _{SM}* to the functional form

which asymptotes to *D _{SM}* as

*T*increases. Here the “best performing RPP model” refers to the RPP model that most accurately reproduced the self-diffusion coefficient computed from 64 particle KS-MD simulations. The parameters

*a*and

*b*are reported in Table III for C at 2.267 g/cm

^{3}and V at 6.11 g/cm

^{3}, and their values vary considerably between both cases emphasizing the need for a comprehensive finite size effect study to produce correction factors for additional elements and conditions. This correction factor allows for the use of the Stanton–Murillo model in regions of previously unknown accuracy. The finite-size corrections along with the corrected Stanton–Murillo model results are shown Fig. 14 with the numerical values given in Table II. Note that for low temperature C at 2.267 g/cm

^{3}, the best performing RPP model was NPA (as reported in Fig. 10 and Table I) explaining why the corrected Stanton–Murillo model tends toward the NPA RPP at low temperatures. For V at 6.11 g/cm

^{3}, the best performing RPP model was the force-matched RPP again explaining the low temperature trend.

Element . | a
. | b
. |
---|---|---|

C (2.267 g/cm^{3}) | 2.198 | −1.032 |

V (6.11 g/cm^{3}) | 0.037 67 | −0.311 2 |

Element . | a
. | b
. |
---|---|---|

C (2.267 g/cm^{3}) | 2.198 | −1.032 |

V (6.11 g/cm^{3}) | 0.037 67 | −0.311 2 |

In an attempt to summarize our work in a single figure, Fig. 15 shows our suggested use cases for all RPPs studied here for two relative self-diffusion accuracies computed from Table I. When points (the average value or its uncertainty) for a given model are within the appropriate tolerance (30% for the top panel and 15% for the bottom panel), we consider the model as being accurate for that temperature and element and is denoted with a colored bar or arrow. We rank the computational expense from lowest to highest as TFY, NPA, force matching, and KS-MD. When a computationally cheaper model is accurate, it replaces the more computationally expensive model in Fig. 15. Based on trends observed in Figs. 6, 11, and 14, we assume that the models remain accurate for higher temperatures and illustrate this by upward pointing colored arrows. For example, consider the case of Fe in the top panel of Fig. 15. The force-matched RPP is accurate to within 30% of the KS-MD result from T = 0.5 eV and up. The NPA model, which is computationally cheaper than the force-matched RPP, becomes accurate (within 30% of KS-MD) at *T *=* *2 eV and up, hence the transition between the force-matched and NPA models. For Al, the NPA RPP is within 15% of KS-MD at all temperatures. However, at *T *=* *15 eV the TFY model becomes accurate therefore replacing the NPA RPP.

## IV. CONCLUSIONS AND OUTLOOK

A systematic study of various RPPs for molecular dynamics simulations of dense plasmas was performed for a wide range of elements vs temperature for solid and half-solid density cases. Of the RPPs studied in this work, RPPs constructed from a NPA approach come closest to accurately reproducing the transport and structural properties predicted by KS-MD. The failures of NPA for metals near *T *=* *0.5 eV are expected: V is a polyvalent metal and s-d hybridization occurs in Au, which is not treated at all in our variant of the NPA model. Thus, it is unclear if inaccuracies in NPA reveal the need for *N*-body interactions or an improved NPA treatment. Moreover, finite-size corrections to KS-MD are seen to be significant; prior work on Si suggests that at least 108 particles are needed to accurately treat elements like C at low temperatures.^{59} Studies on C and Si where there are transient covalent bonding at low temperatures have raised the inadequacy of the PBE XC-functional that has been used here. In Ref. 59, the SCAN functional was used showing remarkable agreement between VASP calculations and NPA results for supercooled high-density Si. This implies that VASP calculations for systems in the low temperature warm dense matter regime become sensitive to the choice of the XC-functional. Similarly, the XC-functionals for transition metals like V, Fe, etc., are known to need Hubbard-type corrections that are not included in our studies. Although this work does not fully resolve these issues, the trends seen for the lowest temperature for C, V, and Au should be examined in detail in future work. Additionally, the NPA model is exceptionally accurate for Al. As Al is a free electron metal, its electronic structure is well described as a Fermi-liquid, the precise physical model in which NPA performs well. In the cases where the electronic structure of the system is not well described as a Fermi-liquid, the performance of the NPA model decreases at low temperature, further emphasizing the need for a comprehensive study over a range of elements and conditions.

As in previous works,^{9,60} the TFY model predicts the least structured *g*(*r*). Notionally, the accuracy of the TFY model appears to follow the machine learning trend of $\u27e8Z\u27e9/Z>0.35$^{61} although it was not possible to use all models in this work at high enough temperatures to be quantitative. In contrast, the NPA model with its improved Kohn–Sham treatment and use of a pseudopotential in (7) eliminates most of these errors except for C and V at *T *=* *0.5 eV, elements for which we would recommend NPA for *T *>* *2 eV. Because we examined seven diverse elements over the warm dense matter regime, the accuracy of NPA (and for moderate temperature, even TFY) suggests that no additional “short-range repulsion”^{9,47–52} is needed beyond (7); as (7) does not contain core–core repulsion, the structure of the interaction is more likely to be effective core–valence repulsion captured by $uei(k)$, as well as structure in $\chi (k)$ beyond $\chi TF(k)$. However, we note (see the Appendix) that in treating weakly ionized systems like warm-dense Ar with a mean ionization of $\u27e8Z\u27e9=0.3$, some 70% of the Ar atoms are neutral, while about 30% of the atoms are singly ionized. Thus, the neutrals interact via a core–core interaction screened by the free electrons. In such cases the use of (7) alone is inadequate. The NPA model treats such a two-component mixture using three pair potentials. In general, core–core interactions are important for weakly ionized atoms with a large core. These core–core interactions can be readily calculated using the core–electron density obtained from the NPA Kohn–Sham calculation.

As expected, the force-matched RPP reproduced the *g*(*r*) computed from KS-MD for all cases. In only one case, again C at 2.267 g/cm^{3} and *T *=* *0.5 eV, the force-matched RPP overestimated the self-diffusion coefficient; this suggests that the spherical pair interaction is not applicable, and non-spherical corrections, which could include three-body contributions, are needed as suggested by the near-perfect agreement of the SNAP and KS-MD microfield of force magnitudes in Fig. 5. However, for all cases considered with *T *>* *1 eV, the *g*(*r*) and self-diffusion coefficient are adequately described by a RPP. With the force-matched-validated NPA interaction, pre-computing the interaction allows for much larger pair-potential simulations.

As fast analytic expressions for transport coefficients are needed for hydrodynamic modeling, we compared our self-diffusion results from all models to the Stanton–Murillo model for both C and V. In both cases, the Stanton–Murillo model was consistent with the TFY model (on which it is based) and both have agreement with force-matched-based results. The error between the Stanton–Murillo model and the force-matched results is $<65%$ below *T *=* *10 eV for V and $<25%$ below *T *=* *5 eV for C, adding confidence to the use of this model in hydrodynamics models above that temperature. For experiments that are rapidly heated above a few eV, little time is spent where the errors are large; because the transport coefficients are numerically very small during this transient heating, negligible transport can occur during that time. For example, note that the V diffusion coefficient varies by a factor of about 30 in the range *T *=* *0.5–100 eV. Conversely, for experiments that dwell at lower temperatures, we provide a RPP-based correction factor to the Stanton–Murillo model with an error of less than 1% for C at *T *=* *0.5 eV and 6% for V at *T *=* *0.5 eV.

Our results suggest several new avenues of investigation. From a data science perspective, larger collections of systematically obtained simulation results would aid in better defining accuracy boundaries. In particular, more elements that produce more material types should be studied. For mixtures, *N*-body potentials could be explored; here, we cast all of the pair potentials as heteronuclear. Additionally, our conclusions are based on studies of the microfield distribution of forces, Einstein frequency, power spectrum, self-diffusion coefficient, and *g*(*r*), which could be extended to include other properties such as viscosities and interdiffusion in mixtures, electrical conductivity, thermal conductivity, and ion-dynamical properties like the speed of sound.^{15} While in this work, we focused primarily on force matching, effective interaction potentials can be obtained through “structure matching.”^{30,62,63} Finally, as very large scale simulations become more common, spatially heterogeneous plasmas can be modeled; much less is known about potentials in such environments, although recent work has explored non-spherical potentials.^{1}

## ACKNOWLEDGMENTS

Sandia National Laboratories is a multimission laboratory managed and operated by 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. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. The authors would like to thank Jeffrey Haack (LANL), Liam Stanton (SJSU), Patrick Knapp (SNL), and Stephanie Hansen (SNL) for insightful discussions, and Josh Townsend (SNL) for his library of VASP post-processing tools. LAMMPS can be accessed at http://lammps.sandia.gov.

## DATA AVAILABILITY

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

### APPENDIX: NPA DETAILS AND EXTENSIONS

In this appendix, we provide a brief background on NPA calculations, and specific extensions needed for the argon, iron and vanadium cases at low temperatures. While we believe the results are formally correct for these cases, we address the nature of the questions being asked and the extensions that are needed.

##### 1. NPA formulation

The NPA model applies primarily to warm-dense fluids with spherical symmetry although the NPA method can be applied equally well to crystals.^{64,65} Importantly, there is no unique NPA model;^{40,42,66,67} here, we describe a specific set of choices^{14,65,68,69} based around a formal statement of the theory.^{69} A key difference between many average-atom models and the NPA model used here is that the free electrons are not confined to the Wigner–Seitz sphere, but move in all of space as approximated by a very large correlation sphere, of radius *R _{c}* which is ten to twenty times the Wigner–Seitz radius.

^{40}

Our NPA model begins with the variational property of the grand potential $\Omega [ne,ni]$ as a functional of the one-body densities $ne=ne(r)$ for electrons, and $ni=ni(r)$ for ions. Only a single nucleus of the material is used and taken as the center of the coordinate system. The other ions (“field ions”) are replaced by their one-body density distribution $\rho (r)$: DFT asserts that the physics is solely given by the one-body distribution; i.e., we do not need two-body, three-body, and such information as they get included via exchange-correlation (XC)-functionals. That is, there is no cluster expansion of the total potential (1) as in (2) of the main text. The terms beyond the pair-interaction are not neglected, but included in the ion–ion XC potential which is not used in VASP-type calculations. Note that this formulation differs from *N*-center codes^{24,70} like the VASP or ABINIT. Moreover, there can be other differences; for example, we have used the finite-*T* electron XC-functional by Perrot and Dharma-wardana,^{71} while the PBE implemented on VASP is a *T *=* *0 XC-functional. The finite-*T* functional used is in good agreement with quantum Monte Carlo XC-data^{68} in the density and temperature regimes of interest.

The artifice of using a nucleus at the origin converts the one body ion density $\rho (r)$ and the electron density *n*(*r*) into effective two body densities in the sense that

The origin need not be at rest; however, most ions are heavy enough that the Born–Oppenheimer approximation is applicable. Here $n\xafi,n\xafe$ are the mean ion density and the mean free electron density, respectively. Bound electrons are assumed to be firmly associated with each ionic nucleus and contained in their “ion cores” of radius *r _{c}* such that

In some cases, e.g., some transition metals, and for continuum resonances etc., this condition for a compact core may not be met, and additional steps are needed. We assume a compact core as a working hypothesis. The DFT variational equations used here are

These directly lead to two coupled Kohn–Sham equations where the unknown quantities are the XC-functional for the electrons, and the ion-correlation functional for the ions.^{72} If the Born–Oppenheimer approximation is imposed, the ion–electron XC-functional may also be neglected. Approximations arise in modeling these functionals and in decoupling the two Kohn–Sham equations^{14,73} to some extent, for easier numerical work. The first equation gives the usual Kohn–Sham equation for electrons moving in the external potential of the ions. This is the only DFT equation used in *N*-center codes in which ions define a periodic structure evolved by MD, followed by a Kohn–Sham solution at each step. In contrast, NPA employs the one-body ion density $ni(r)$; it was shown in Ref. 69 that the ionic DFT equation can be identified as a Boltzmann-like distribution of field ions around the central ion, distributed according to the “potential of mean force” well known in the theory of fluids. In such a formulation, the ion–ion correlation functional $Fxcii$ was identified to be the sum of hypernetted-chain diagrams plus the bridge diagrams as an exact result formally although the bridge diagrams cannot be evaluated exactly.

The mean electron density $n\xafe$ can also be specified as the number of free electrons per ion, viz., the mean ionization state $\u27e8Z\u27e9$. Although the material density $n\xafi$ is specified, the mean free electron density $n\xafe$ is unknown at any given temperature, as it depends on the ionization balance which is controlled by the free energy minimization given in (A3). Hence, a trial value for $n\xafe$ (i.e., equivalently, a trial value for $\u27e8Z\u27e9$) is assumed and the thermodynamically consistent $ni(r)$ is determined. This is repeated until the target mean ion density $n\xafi$ is obtained.

This means that the Kohn–Sham equation has to be solved for a single electron moving in the field of the central ion; its ion distribution $\rho \xafg(r)$ is modified at each iteration with modification of the trial $Z\xaf$ until the target material density is found. However, it was noticed very early^{14,65} that the Kohn–Sham solution was quite insensitive to the details of the *g*(*r*) and hence a simplification was possible. The simplification was to replace the trial *g*(*r*) at the trial $Z\xaf$ by a cavity-like distribution

Here the *a _{i}* is the trial value of

*a*, based on the trial $n\xafe$. Hence, adjusting the $gcav(r)$ at each iteration requires only adjusting the trail

_{i}*a*to achieve self-consistency. The self-consistency in the ion distribution is rigidly controlled by the Friedel sum rule for the phase shifts of the Kohn–Sham-electrons.

_{i}^{69}This ensures that $\rho \xaf=n\xaf/Z$. Thus, a valuable result of the calculation using the cavity model of

*g*(

*r*) is the self-consistent value of the mean ionization state

*Z*, which is both an atomic quantity and a thermodynamic quantity.

Here we note crucial simplifications used in implementing the NPA. Given that the electron distribution *n*(*r*) obtained self-consistently can be written as a bound-electron term and a free-electron term because of the condition specified by (A2), we have

The core–electron density (made up of “bound electrons”) is denoted by *c*(*r*). The free electron density $nf(r)$ is the response of an electron fluid containing a cavity that mimics $ni(r)$. It contributes to the potential on the electrons. The response of a uniform electron gas to the central ion can be obtained by subtracting the effect of the cavity using the known static interacting linear response function $\chi (k,n\xafe,T)$ of the electron fluid. That is, from now on we take it that the charge density $nf(r)$ and the charge pileup $\Delta nf(r)$ are both corrected for the presence of the cavity, but we use the same symbols.

##### 2. Pair-potentials for NPA mixtures: Ionic contributions

With the basic NPA formulated, we turn to the construction of pair potentials, with a focus on the more challenging cases we consider in the main text. It is very common to treat WDM and liquid metals through purely ionic interactions, which are adequately evaluated in second-order perturbation theory unless the free electron density and the temperature ($T/EF$) are very low. Such interactions, which generalizes (7) to mixtures, so are written in *k*-space as

In the NPA theory for mixtures, *Z _{a}* and

*Z*are integers, while in the simple (average-atom) NPA, the $Zs\xaf=\Sigma sxsZs$ is used. The electron density pileup is calculated using the linear-response property of the pseudopotential

_{b}Here, since $nf(k)$ has been calculated via Kohn–Sham, it has all the non-linear effects included by the construction of $Uei(k)$. The extent of the validity of such a quasi-linear pseudopotential is discussed in Ref. 13. Unlike in the average atom NPA, the *U ^{ei}* used in mixture-theory is the pseudoptential of the ion with the appropriate

*integer ionization*. The interacting electron gas response function used in these calculations includes a local-field factor chosen to satisfy the finite temperature electron-gas compressibility sum rule, and is given explicitly by

Here *χ*_{0} is the finite-*T* Lindhard function, *V _{k}* is the bare Coulomb potential, and

*G*is a local-field correction (LFC). The finite-

_{k}*T*compressibility sum rule for electrons is satisfied since

*κ*

_{0}and

*κ*are the non-interacting and interacting electron compressibilities, respectively, with

*κ*matched to the $Fxc(T)$ used in the Kohn–Sham calculation. In (A12), $kTF$ appearing in the LFC is the Thomas–Fermi wave vector. We use a

*G*evaluated at $k\u21920$ for all

_{k}*k*instead of the more general

*k*-dependent form [e.g., Eq. (50) in Ref. 71] since the

*k*-dispersion in

*G*has negligible effect for the WDMs of this study.

_{k}##### 3. Pair-potentials for NPA mixtures: Core interactions

We now consider extensions to this formulation that includes the effect of core electrons (see Appendix B of Ref. 14). Core–core interactions are important for atoms like argon, sodium, potassium, gold, etc., with large cores and zero or low $\u27e8Z\u27e9$. Here we use a simplified approach in discussing the case of low-*T* argon rather than the more exact approach given in Ref. 14. The total pair-interaction *ψ _{ab}* is of the form

The first term, $Uc=U(ca,cb)$ is the interaction between the two atomic cores, and it is the only term found in a neutral gas of pure argon atoms. The second term is the interaction of the core electrons of the neutral atom *a* with the pseudo-potential of the ion *b* with integer charge *Z _{b}*, while the indices are interchanged in the third term. We need to evaluate the first two terms involving core electrons, while the last term is given by (7) of the main text. How to evaluate these in the NPA at any temperature and density is given in Appendix B of Ref. 14 to the same level of approximation as (A6), i.e., to second-order in perturbation theory in the screened interactions. It is found that such evaluations for argon give results close to parameterized potentials similar to the Lennard–Jones (LJ) or more sophisticated potentials, but inclusive of a screening correction. Thus, at the LJ-level of approximation, $Uc(k)$ for two neutrals immersed in the electron gas is approximately given by $ULJ(k){1+Vk\chi (k)}$. For two charged ions, a correction factor of ${(Zn\u2212Zint)/Zn}2$, where

*Z*is the nuclear charge, and

_{n}*Z*is the integer-ionization of the ion is included because the ion cores have less electrons than the neutral cores.

_{int}##### 4. Argon in low-T WDM states

The NPA calculation at density $\rho \u2243$ 1.4 g/cm^{3} and temperature *T *=* *2 eV yields a mean charge of $\u27e8Z\u27e9\u22430.3$; however, this implies a mixture where 30% of the argon atoms are in the Ar^{+} state, while 70% of the atoms are neutral atoms. We ignore the Ar^{2+} ionization state as the second ionization energy is about 27 eV (ignoring plasma corrections). So, while the NPA tries to assign a mean ionization state, this argon system is better treated as a two component mixture with *x _{a}*,

*x*, where $a=Ar,b=Ar+$.

_{b}Thus, we need the three pair-potentials *U _{ab}*, i.e., for Ar–Ar, Ar–Ar

^{+}, and Ar

^{+}–Ar

^{+}interactions. Here the atomic species are immersed in the electron gas resulting from the ionization. These pair-interactions can be rigorously calculated from first-principles using the atomic and ionic electron densities obtained in the NPA, as discussed in Appendix B of Ref. 14. Here we follow a more simplified calculation of the pair potentials exploiting known models for argon, suitable for this mixture with 70% neutral argon.

The Ar^{+} ion interacts with the neutral Ar atom by polarizing the core–electron distribution. This distortion is usually described in terms of the dipole polarizability *α* and quadrupole polarizability *β* of the Ar atom, viz.,

This interaction will be screened by the polarization processes of the background electron gas. Instead of using the truncated multipole expansion, an improved calculation may be done as in Appendix B of Ref. 14. The ion–atom contribution from (A14) contains only the first term, as atom *a* is neutral and has no free electrons. In the regime *T *=* *2 eV for Ar, we have simply approximated the core–core interaction of the atom–ion interaction by the mean value of atom–atom and ion–ion core–core interaction.

In Fig. 16, we display the pair potentials for Ar–Ar, Ar^{+}–Ar, and Ar^{+}–Ar^{+} at *T *=* *2 eV at a density of 1.395 g/cm^{3}. The compositions 0.7 for Ar, and 0.3 for Ar^{+} are obtained from the simple NPA calculation. We have not attempted to optimize these composition fractions using steps indicated in Ref. 14.

##### 5. Iron and Vanadium in low-T WDM states

The NPA uses an “isotropic” atomic model even for iron, vanadium, and other transition metals, and ignores the *s*-*d* hybridization energy *E _{hyb}* that re-arranges the electron distribution among the

*d*and

*s*electron states near the Fermi energy. The need for

*s*-

*d*hybridization is best seen by looking at the low temperature band structure of such a transition metal. The unhybridized bands of an “isotropic model” for vanadium are such that the free-electron like band crosses the

*d*-bands, and the

*s*-

*d*interaction redistributes the electrons. Furthermore, instead of the

*s*-wave local pseudopotential (as employed here), an angular-momentum dependent form is more appropriate at low-

*T*. Hence, the calculation of the mean ionization has to be accordingly modified. However, the simpler picture reemerges when the temperature exceeds the

*s*-

*d*hybridization energy.

Furthermore, charge polarization fluctuations of the *d*-electrons can couple with those of the electron gas (e.g., as discussed by Maggs and Ashcroft^{74}). Such effects, as well as “on-site Hubbard U” effects need to be included in the electron XC-functionals to properly treat transition metals. The usual XC-functionals made available with standard codes do not include these effects. However, at sufficiently high temperatures, the ionization becomes high enough to screen these effects and the theory simplifies. The pair-potentials for iron provided by the isotropic model with no *s*-*d* hybridization and other *d*-band effects lead to cluster formation at low-*T*. Here, unlike in liquid carbon or silicon, the bonding is not transient. Hence low-*T* MD simulations will show no movement of the ions and no diffusion. In the case of silicon and carbon, both of which are known to show varying degrees of transient bonding behavior, NPA was shown to produce an accurate description of the high density liquid phase, with structure factors and PDFS obtained from the NPA-HNC procedure agreeing closely^{75} with those from 218-atom simulations with DFT.