Efficacy of the Radial Pair Potential Approximation for Molecular Dynamics Simulations of Dense Plasmas

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 $\sim20\%$.


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. 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 * staneklu@msu.edu † rclay@sandia.gov ‡ murillom@msu.edu In this work, we carry out KS-MD simulations for a range of elements, temperature, and density, allowing for a systematic comparison of three RPP models. While multiple RPP models can be selected, 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 [2][3][4][5], and the optimal force-matched RPP that is constructed directly from KS-MD simulation data.
With this choice of RPP models, we constrain the functional form of the RPP required to accurately describe microscopic interactions for the systems studied here. For example, if the forced-matched RPP cannot accurately reproduce the metrics described herein, it is very likely that no RPP can describe the system and a KS-MD simulation, or perhaps an interaction potential beyond a RPP (e.g. three-body interaction potentials) is needed. Additionally, if NPA results agree with KS-MD, we are able avoid the force-matching process which requires data from KS-MD. Furthermore, if a Yukawa model agrees with KS-MD, we can neglect including a pseudopotential in the construction of the RPP, reducing the interaction to a simple, analytic form.
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 g(r), Einstein frequency, power spectrum, and self-diffusion coefficient 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 emphasises 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, n i , 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 ρ i , or in some cases half of that, and in thermodynamic equilibrium at temperature 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 U tot = U N (r 1 , r 2 , . . . , r N ). (1) 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 N o × N o 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 onebody problem, matrix diagonalization incurs a cost of O(N 3 o ), 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 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, n i , T }. The second approach to calculating the po-tential energy surface is to use a cluster-type expansion, which takes the form (2) 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 [6][7][8]. 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 [9][10][11], 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 f i (T ) = (1 + e β(Ei−µ) ) −1 , and the Kohn-Sham-Mermin orbitals φ 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) [12][13][14][15]. The finite temperature electronic structure was treated with Mermin free-energy functional, and we used the Perdew-Burke-Ernzerhof functional for the exchange correlation energy [16]. To improve computational efficiency, we eliminated the chemically inactive core electrons with projector augmented-wave [17] pseudopotential, the details of which are given in the supplemental information [18]. 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 [19] 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 was 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 an approximate functional form for (2) and defining a suitable cost-function based on the difference between the Kohn-Sham and approximate model forces, we use the force matching procedure [20][21][22][23][24] to generate the optimal RPP model based on the KS-MD force data. By choosing an optimizable spline 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 3-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 [25,26]. 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 [27] written in terms of the mean ionization state Z , the bare Coulomb potential u C (k) = 4πe 2 /k 2 , the electron-ion pseudopotential u ei (k) and the susceptibility χ(k).
In practice, pair interactions are constructed using the 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 [28]; such calculations can also be relativistic [29,30]. From the AA, a pseudopotential u ei (k) and an accurate free/valence electron response function χ(k) are constructed and (6) is formed. This approach has three strengths: (1) typical AA models are not limited to low temperatures, (2) the interaction (6) 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 a i = (3/4πn i ) 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 (6) one approximates the pseudopotential as u ei (k) ≈ −4π Z e 2 /k 2 , where the mean ionization state Z results from a AA calculation [31], and χ(k) in its long-wavelength (Thomas-Fermi) limit χ T F (k); this is known as the "Yukawa" interaction [32,33]. 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 and the Fermi energy E F = 2 (3π 2 n e ) 2/3 /2m e . 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 u ei (k) and χ(k). The TFY model can be extended by including the gradient corrections to χ T F (k), but otherwise retaining the other approximations. This improvement yields the Stanton-Murillo potential [33]; the gradient correction to χ T F (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 (6). 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 [34]. 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 exchangecorrelation potential; the susceptibility is Lindhard with local field corrections [3]. Note that the electron-ion pseudopotential u ei (k) introduces additional oscillations on length scales different from χ(k), although the Friedel oscillations in χ(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 > E F ; 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 (6) 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/r 4 [35][36][37][38][39][40][41]; 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/r 4 is only valid over a very small range of r values; importantly, the NPA interaction shows that a increases with r, 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 [2]. Worse, the short-range repulsion model potentially gets an accurate answer for the wrong reason, as we explore in Fig. 3.
Because the form (6) 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 (6) naturally contains both the "crowding" and "attraction" behaviors as special cases. Fig. 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 attractive have attractive 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 Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [42]. 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 3 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.

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 an approximate model (RPP or many-body potential) for M particle coordinate configurations. For this, we compute the mean-absolute force error are the α-th force components (x, y, or z) on the i-th atom in particle coordinate configuration number m for the approximate 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 distribution 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 predicts more frequent small force magnitudes relative to force magnitudes computed from the KS-MD force data. In contrast, for C, V, and Au at T = 0.5 eV the NPA RPP tends to predict more large force magnitudes relative to the KS-MD force data. These trends can be connected back to (6) where the choice of Z , u ei (k), and χ(k) all contribute to the construction of a RPP model and hence the force magnitudes. More work needs to be done FIG. 4. The RPP models normalized by temperature versus distance for C, Al, V, and Au. Top row, T = 0.5 eV: the representation of the RPP is element dependent with strong agreement for aluminum. Bottom row, T = 5 eV: The agreement between models improves significantly. The differences in the representation can be connected back to (6) where the treatment of the mean ionization, electron-ion pseudo potential, and susceptibility define the RPP.  Table I. In general, when the microfields are similar to that of KS-MD, the agreement between the self-diffusion coefficient increases. To assess the importance of 3-body or higher interactions, SNAP results are reported for C and V at both T = 0.5 and 5 eV.
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 (9) to normalize the results across the different elements, temperatures, and densities studied here. Such a scale factor is the "mean absolute force" defined as Using (9) and (10) we define the relative force error 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, (11) provides a temperature independent metric as measured with respect to a KS-MD force data "baseline". Intuitively, when (11) evaluates to 1, the mean absolute error is the same order of magnitude as the mean absolute force and when (11) is 0, the approximate model is exactly reproducing the per-component KS-MD force data. Fig. 6 displays (11) 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 decrease towards 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 3-body and higher interactions by including relative force errors using a SNAP model. We observe that the relative force error decreases significantly relative to the force-matched RPP. For C, the relative force error drops from roughly 0.6 us-ing 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 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, which reproduce the KS-MD reference forces to within a comparable amount relative to other elements, produce a self-diffusion coefficient that differs from the KS-MD result by many factors. 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 force error is high. The imperfect mapping of relative self diffusion error versus 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 measures the probability of finding another ion at a radial distance r around a central ion. It has been shown that in general, there always exists a RPP that can reproduce g(r) from a Nbody simulation [43] and the force-matching procedure provides a avenue for obtaining this RPP. Fig. 8 compares g(r) computed from MD simulation 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.
To understand how underestimating the g(r) peak height impacts the self-diffusion coefficient, we use the Einstein frequency which is obtained through a short time expansion of the normalized velocity autocorrelation function where Ω 0 is the Einstein frequency where m i is the ion mass in grams. The Einstein frequency gives insight on 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 (13) is shown. For the TFY model, the integrand is always smaller than those predicted by forcematched 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 3 dimensions) by substituting, (12) into (14). 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
One approach to compute the self-diffusion coefficient is via the the slope of the mean-squared displacement from the Einstein relation where · is an ensemble average over particles and time.
Relative to (14) we have found that (15) converges faster for small numbers of particles and short simulation lengths, and is less sensitive to finite-size effects. Therefore, all self-diffusion coefficients reported in this work have been calculated from a linear fit to the meansquared displacement, |r(t) − r(0)| 2 . 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 latetime 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 new sub-sample fits yields self-diffusion coefficient and the uncertainty in the fit respectively and are reported in Table I. By computing the relative self-diffusion coefficients reported in Table I to determine a rule where NPA is accurate relative to KS-MD and similarly where TFY is accurate relative to NPA. Doing so, we are informed of when certain models may be accurate and when others are not. 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 selfdiffusion 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   Table I as the self-diffusion coefficient decrease as the integral of the Einstein frequency increases. This allows for a "by eye" comparison of the self-diffusion coefficient from different RPP models.  Table I. For all cases all models predict values that have roughly the same order of magnitude. The only case where the force-matched RPP fails to reproduce the KS-MD self-diffusion coefficient is for C at 2.267 g/cm 3 and T = 0.5 eV. General trends in the data include the over prediction in the self-diffusion coefficient from the TFY RPP, consistent with the Einstein frequency in Fig. 9. This over prediction is reduced greatly when using a NPA RPP and is generally in good agreement with the force-matched RPP and KS-MD data.
NPA is adequate. The bottom panel shows a comparison with TFY, which has the simplest u ei (k) and χ(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 compareZ(ν) calculated using TFY, forcematched, 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 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 highfrequency 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 is agreement comes at the expense of sacrificing the accuracy of static properties. Lastly, the TFY RPP exhibits roughly the same trends across all elements and temperaturesoverestimation 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,6]. 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 [44] where D ∞ is the self-diffusion in the thermodynamic limit, D N is the self-diffusion coefficient computed from a system of finite number of particles N , and ξ = 2.837297 for cubic simulation boxes 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 ∞ . Results from this procedure are shown in Fig. 13 where D ∞ is determined via linear extrapolation to 1/L = 0. By finding the percent difference in D ∞ 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 ∼ 20%. 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 1 eV).
Finite-size corrections allow for a direct comparison to analytic transport theories, namely the Stanton-Murillo model [45]. The Stanton-Murillo model, provides a closed form solution for ionic self-diffusion by using an effec-tive 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 [45], more work is needed to quantify these trends for other transport properties.
With the converged self-diffusion data, we generate an effective interaction correction to original Stanton-Murillo model.
The effective interaction corrected Stanton-Murillo model is where α(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 towards 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.
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 FIG. 12. The normalized power spectrum for C, Al, V, and Au. For C at T = 0.5 eV, the single particle dynamics are poorly described by the TFY and force-matched models but more accurately described with the NPA model. As the temperature increases from T = 0.5 to 5 eV, all models more accurately reproduce small and high frequency dynamics with the most notable improvement for C. TABLE II. Self-diffusion coefficients in the thermodynamic limit. Both elements are at solid density (2.267 g/cm 3 for C, and 6.11 g/cm 3 for V). 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 a models remain accurate for higher temperatures and is illustrated by the upward pointing colored arrows. 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.
FIG. 13. Finite-size effect study for V at 6.11 g/cm 3 and T = 2 eV. Identical MD simulations were carried out with increasing particle number. Extrapolating with a linear fit (grey dashed line) to 1/L = 0 approximates the thermodynamic limit, correcting the values in Table I.   TABLE III. Coefficients a, and b for the effective interaction correction (19). Note that the values of a and b vary considerably for each element.

IV. CONCLUSIONS AND OUTLOOK
A systematic study of various RPPs for molecular dynamic simulations of dense plasmas was performed for a wide range of elements versus 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 are not treated well 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 Silicon suggests that at least 108 particles are needed to accurately treat elements like C at low temperatures [46]. 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.
As in previous work [35,47], the TFY model tends to underestimate the peak (relative to KS-MD) of g(r) at FIG. 14. Self-diffusion coefficient versus temperature in the thermodynamic limit. The points displayed here are taken from Table II. (a) Self-diffusion coefficient for C at 2.267 g/cm 3 .(b) Self-diffusion coefficient for V at 6.11 g/cm 3 . The Stanton-Murillo model (denoted SM) fails for low temperature C. For V, the Stanton-Murillo model shows excellent agreement with the force-matched RPP even at low temperatures. The validity of the Stanton-Murillo model is extended to low temperatures with an effective interaction correction (denoted CSM). low temperatures, with a corresponding error in the self diffusion coefficient. Notionally the accuracy of the TFY model appears to follow the machine learning trend of Z /Z > 0.35 [48], although it was not possible to use all models at high enough temperatures to be quantitative. In contrast, the NPA model with its improved Kohn-Sham treatment and use of a pseudopotential in (6) 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" [35][36][37][38][39][40][41] is needed beyond (6); as (6) does not contain corecore repulsion, the structure of the interaction is more likely to be effective core-valence repulsion captured by u ei (k), as well as structure in χ(k) beyond χ T F (k). In fact, our results suggest that core-core repulsion is quite small as ions do not approach each other on core length scales because of Coulomb repulsion; at higher temperatures, where ions could have closer encounters, the core is smaller; thus, the first term in (6) is purely Coulombic.
FIG. 15. Suggested use cases for each RPP model based on the relative self-diffusion coefficient error (between RPP-MD and KS-MD) and cheapest computation cost. The top and bottom panels correspond to a 30% and 15% relative error respectively. The elements denoted with a subscript of "1/2" corresponds to half solid density (V at 3.055 g/cm 3 and Fe at 3.937 g/cm 3 ). The colored bars indicate the computationally cheapest RPP that generates a self-diffusion coefficient to within the specified error tolerance available for that system based on Table I. The empty space under each bar indicates regions where no KS-MD data was collected so no assessment on a RPPs accuracy can be made.
As expected, the force-matched potential 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 forcematched RPP overestimated the self diffusion coefficient; this suggests that the spherical pair interaction isn't 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, precomputing the interaction allows for much larger pairpotential 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 Stanton-Murillo error relative to force-matched 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 to 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 g(r) and self-diffusion, which could be extended to include other properties such as viscosities and interdiffusion in mixtures. Finally, as very large scale simulations become more common, spatially heterogeneous plasmas can be modeled; much less is known about potentials of any order in such environments, although recent work has explored non-spherical potentials [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 [49,50]. Importantly, there is no unique NPA model [31,[51][52][53]; here, we describe a specific set of choices [50,[54][55][56] based around a formal statement of the theory [56]. A key difference between many average-atom models and the NPA 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 [31].
Our NPA model begins with the variational property of the grand potential Ω[n e , n i ]) as a functional of the onebody densities n(r) for electrons, and ρ(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 ρ(r): DFT asserts that the physics is solely given by the one-body distribution; i.e., we do not need twobody, three-body, and such information as they get included via exchange-correlation (XC)-functionals. Note that this formulation differs from N -center codes [57,58] 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 [59], while the PBE implemented on VASP is a T = 0 XCfunctional. The finite-T functional used is in good agreement with quantum Monte-Carlo XC-data [54] in the density and temperature regimes of interest.
The artifice of using a nucleus at the origin converts the one body ion density ρ(r) and the electron density n(r) into effective two body densities in the sense that n i (r) =n i g ii (r), n e (r) =n e g ei (r) (A1) The origin need not be at rest; however, most ions are heavy enough that the Born-Oppenheimer approximation is applicable. Heren i ,n e 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. These directly leads to two coupled Kohn-Sham equations where the unknown quantities are the XCfunctional for the electrons, and the ion-correlation functional for the ions [60]. 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 [55,61] 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 n i (r); it was shown in [56] 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 F ii xc 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 densityn e can also be specified as the number of free electrons per ion, viz., the mean ionization state Z . Although the material densityn i is specified, the mean free electron densityn e is unknown at any given temperature, as it depends on the ionization balance which is controlled by the free energy minimization A3. Hence, a trial value forn e (i.e, equivalently, a trial value for Z ) is assumed and the thermodynamically consistent n i (r) is determined. This is repeated until the target mean ion densityn i 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ρg(r) is modified at each iteration with modification of the trialZ until the target material density is found. However, it was noticed very early [50,55] 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 trialZ by a cavity-like distribution: Here the a i is the trial value of a i , based on the trialn e . Hence, adjusting the g cav (r) at each iteration requires only adjusting the trail a i 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 [56]. This ensures thatρ =n/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 boundelectron term and a free-electron term because of the condition specified by Eq. A2, we have: n e (r) = c(r) + n f (r), ∆n e (r) = n f (r) +n e (A5) The core-electron density (made up of "bound electrons") is denoted by c(r). The free electron density n f (r) is the response of an electron fluid containing a cavity that mimics n i (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 χ(k,n e , T ) of the electron fluid. That is, from now on we take it that the charge density n f (r) and the charge pile up ∆n f (r) are both corrected for the presence of the cavity, but we use the same symbols.

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 /E F ) are very low. Such interactions, which generalizes (6) to mixtures, so are written in k-space as In the NPA theory for mixtures, Z a and Z b are integers, while in the simple (average-atom) NPA, thē Z s = Σ s x s Z s is used. The density fluctuations use the linear-response property of the pseudopotential Here, since n f (k) has been calculated via Kohn-Sham, it has all the non-linear effects included by the construction of U ei (k). The extent of the validity of such a quasi-linear pseudopotential is discussed in Ref. [62]. 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 k is a local-field correction (LFC). The finite-T compressibility sum rule for electrons is satisfied since κ 0 and κ are the non-interacting and interacting electron compressibilities respectively, with κ matched to the F xc (T ) used in the Kohn-Sham calculation. In Eq. A12, k TF appearing in the LFC is the Thomas-Fermi wave vector. We use a G k evaluated at k → 0 for all k instead of the more general k-dependent form (e.g., Eq. 50 in Ref. [59]) since the k-dispersion in G k has negligible effect for the WDMs of this study.

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. [55]). Core-core interactions are important for atoms like argon, sodium, potassium, gold, etc., with large cores and zero or low Z . Here we use a simplified approach in discussing the case of low-T argon rather than the more exact approach given in Ref. [55]. The total pairinteraction ψ ab is of the form: The first term, U c = U (c a , c b ) 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. The last term is the corresponding DFT corrections to the exchange-correlation energy and to the interacting kinetic energy. How to evaluate these in the NPA at any temperature and density is given in the Appendix B of Ref. [55] to the same level of approximation as Eq. 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 parametrized potentials similar to the Lennard-Jones(LJ) or more sophisticated potentials, but inclusive of a screening correction. Thus, at the LJ-level of approximation, U c (k) for two neutrals immersed in the electron gas is approximately given by U LJ (k){1 + V k χ(k)}. For two charged ions, a correction factor of {(Z n − Z int )/Z n } 2 , where Z n is the nuclear charge, and Z int is the integer-ionization of the ion is included because the ion cores have less electrons than the neutral cores.

Argon in low-T WDM states
The NPA calculation at density ρ 1.4 g/cm 3 and temperature T = 2eV yields a mean charge of Z 0.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 b , where a = Ar, b = Ar + .
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 reference [55]. 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.
(A14) U (c a , f b )(r) = C 4 /r 4 + C 6 /r 6 + ..., (A15) C 4 = α(Z int ) 2 /2, C 6 = β(Z int ) 2 /2. (A16) 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. [55]. The ionatom contribution from Eq. A14 contains only the last term, as the atom b is neutral and U ei b (k) = 0. 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 Figure A 4 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. [55].

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 re-emerges when the temperature exceeds the s − d hybridization energy.
Furthermore, charge polarization fluctuations of the delectrons can couple with those of the electron gas (e.g., as discussed by Maggs and Ashcroft [63]). 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 [64] with those from 218-atom simulations with DFT.