Simple pair-potentials and pseudo-potentials for warm-dense matter and general applications
aa r X i v : . [ c ond - m a t . m t r l - s c i ] F e b Simple pair-potentials and pseudo-potentials for warm-dense matter applications.
M.W.C. Dharma-wardana ∗ National Research Council of Canada, Ottawa, Canada, K1A 0R6 (Dated: February 23, 2021)We present methods for generating computationally simple parameter-free pair potentials usefulfor solids, liquids and plasma at arbitrary temperatures. They successfully treat warm-dense matter(WDM) systems like carbon or silicon with complex tetrahedral or other structural bonding features.Density functional theory asserts that only one-body electron densities, and one-body ion densitiesare needed for a complete description of electron-ion systems. DFT is used here to reduce both the electron many-body problem and the ion many-body problem to an exact one-body problem,namely that of the neutral pseudoatom (NPA). We compare the Stillinger-Weber (SW) class ofmulti-center potentials, and the embedded-atom approaches, with the NPA approach to show thatmany-ion effects are systematically included in this one-center method via one-body exchange-correlation functionals. This computationally highly efficient “single-center” DFT-NPA approachis contrasted with the usual N -center DFT calculations that are coupled with molecular dynamics(MD) simulations to equilibriate the ion distribution. Comparisons are given with the pair-potentialparts of the SW, ‘glue’ models, and the corresponding NPA pair-potentials to elucidate how the NPApotentials capture many-center effects using single-center one-body densities. PACS numbers: 52.25.Jm,52.70.La,71.15.Mb,52.27.Gr
I. INTRODUCTION
Condensed-matter systems at temperatures T ≥ per se . Furthermore, ‘warm’ systems are those away fromthe simplifications available for the T → T → ∞ limits. Hence a truly general finite- T many-body the-ory of interacting electrons and ions in arbitrary statesof matter is needed. Density functional theory (DFT)is a favoured approach since it reduces such problems toa theory of mere one-body densities, requiring no multi-atom approaches, at least in principle. Here we exploitDFT for simplifying the general electron-ion problem tothe fullest [1, 3].While typical WDM systems are in the domain ofplasma physics, laser-matter interactions, high pressurephysics, geophysics, astrophysics etc., many problems innano-structure physics also fall into WDM physics. Anelectron layer or a hole layer in a quantum well containsparticles of modified effective mass ( m ∗ ) that may be atenth of the free electron mass. The corresponding Fermienergies and particle densities are such that even at roomtemperature, the electrons may be partially degeneratewhile the holes are classical particles [4]. The quantumwells or nanostructures may have metal-oxide layers andinhomogenieties that need a unified many-body first prin-ciples theory [5].The need for practical computations of energetics ofdefects and dislocations in metallic structures led to semi-empirical theories usually known as embedded-atom or ∗ Email address: [email protected] effective medium theories [9]. The potential V ( ~R ) felt byan ion at ~R is modeled using two-body, three-body andpossibly higher terms in the energy, and used to moti-vate functional forms for numerically fitted representa-tions of structure dependent energies. A parallel devel-opment in semiconductor physics modified the usual two-body potentials of ‘covalent bonds’ (e.g., of the Lenard-Jones type) with ‘structure dependent’ three-body terms,viz., as in the well-known Stillinger-Weber [10] or Ter-soff models [11] for C, Si and other tetrahedral solidsand fluids. These have been extended to include em-pirical ‘bond-orders’, angular and ‘torsional’ correla-tions [12], producing very complex ‘potentials’. Highlyparameter-loaded potentials fitted to empirical data aswell as large DFT data bases [13] have emerged. Manyof these multi-center multi-parameter potentials havebeen incorporated in codes like LAAMPS (Large-scaleAtomic/Molecular Massively Parallel Simulator) [14] forregular use. In these descriptions, the electron systemis fully integrated out and only a classical simulationbased on the multi center potentials is used. Hencethese multi-center effective-medium models lack associ-ated atomic pseudopotentials that may be used to delin-eate their building blocks, or for use in transport stud-ies. In contrast, the NPA potentials that we discuss hererelate directly to underlying pseudopotentials, responsefunctions, and XC-functionals as the NPA is based onfirst-principles physics. It also include the multi-center‘effective medium’ effects via ion-ion XC-functionals gen-erated in situ .Attempts to extend these multi-center potentials, anduse N -center DFT for T > T , with a one-body electron density n ( ~r ), and a one-body ion density ρ ( ~r ), standard DFTasserts that all the thermodynamic properties and lin-ear transport properties are a functional of just the one-body n ( r ) and ρ ( r ), while the many-body interactions arealso functionals of just these one body densities. Hencemulti-center approaches should not be necessary if DFTis rigorously applied. The many-body effects are rele-gated to the exchange-correlation (XC) functionals. Thefunctionals are non-local but remain one-body densityfunctionals. At finite- T we use the free-energy function-als F xcss ′ [ n s , n s ′ ], where s or s ′ = e, i is the species index,with n e = n ( r ) , n i = ρ ( r ) [1]. The finite- T free-energyfunctionals reduce to the usual ground-state energy func-tionals at T = 0.However, common computer implementations of DFTin large codes like ABINIT [16] or VASP [17] explicitlyrequire the locations ~R I , I = 1 , N , of the nuclei of N ions and only the electron problem is reduced by DFT.The density functional theory of the one-body ion den-sity ρ ( r ) is not invoked in these codes. Only an electronexchange-correlation functional F xcee is invoked within aKohn-Sham (KS) calculation for specific ionic configu-ration ~R I using the Born-Oppenheimer approximation.The ions merely provide an ‘external potential’. Classicalmolecular dynamics is used to equilibrate the ion distri-bution P δ ( ~r − ~R I ). Such an approach produces a many-centered potential energy function V ( ~R I ) , I = 1 , · · · , N .So, single-ion properties, pair-potentials between ionsetc., are not directly available by this N -center DFT, asthe calculation is for a solid-state cluster of N -nuclei andelectrons moving on a complex potential-energy surface.Extraction of such ‘atomic’ properties invokes an ex-pansion of, say, the potential felt by an ion at ~R in termstwo-body, three-body and higher interactions. A trunca-tion in finite order and fitting a required property, e.g.,an effective potential, to a large parameter set is used.A kinetic energy functional [18] may be used to simplifythe DFT computations if some accuracy could sacrificed,but ‘one-atom’ properties or pair-potentials are not di-rectly available from such N -ion calculations. Usually N ranges from 100-500 particles unless crystal symmetry orsome such property can be invoked.The advantage of the NPA form of DFT is that single-atom properties, linearized atomic pseudopotentials, aswell as pair-potentials, and the corresponding pair distri-bution functions (PDFs) are directly available and maybe used as inputs to both quantum simulations (using thepseudopotentials) as well as classical simulations usingthe pair-potentials. That these pair-potentials do NOTneed multi-center corrections is not adequately appreci- ated in the community of users. Typical examples ofsuch advantages of the NPA are in providing the meanionization ¯ Z , ionization fractions x s for a mixture of ionsof different ionization states Z s + , providing an electron-ion interaction U ei ( k ), providing PDFs needed for lineartransport calculations, for temperature relaxation stud-ies, providing the ion-feature in X-ray Thomson scat-tering (XRTS), while being able to use Fermi-liquid ap-proaches when needed [20].The usual approach of explicit manipulation of atomicpositions used in large codes [16, 17] is quite appropri-ate for solid-state physics applications based on the unitcell, or for the quantum chemistry of a finite number ofatoms usually at T = 0. They provide microscopic de-tails of ‘bonding’ between atoms not available from theNPA. The NPA gives only a thermodynamic average ofthe system, via the thermodynamic one-body densities n ( r ) , ρ ( r ). In constructing the NPA that reduces the N center problem to a single center, we explicitly identifyonly one nucleus and then consider the one-body distri-butions n s ( r ) , s = e, i around it rather than their explicitpositions, and exploit the spherical symmetry of liquidsand plasmas, or the crystal symmetry of solids, to pro-vide a computationally very efficient and yet completelyrigorous scheme. The generalization to mixtures contain-ing many species increases the number of XC-functionalsneeded, and proceeds as in Ref. [2]. The NPA theory ofmixtures will be illustrated by some explicit examples.Use of DFT for the ion distribution avoids the adhoc introduction of multi-atom effects found in semi-empirical potentials like the Stillinger-Weber potentials,effective medium models, or in the Finnis-Sinclair poten-tials [19] used in metal physics. We have shown thata systematic procedure based on the Ornstein-Zernikeequation exists for including the three-body and higherterms into the ionic XC-functional F xcii . Since ions areclassical particles in most cases, the ‘exchange’ contentof F xcii is negligible.Unlike fitted potentials, the NPA method, being a firstprinciples approach, can be used for unusual states ofmatter. If the electron subsystem is at a temperature T e ,and the ions subsystem is at a temperature T i , we haveone body densities n e ( r, T e ) and ρ ( r, T i ) defining a two-temperature WDM system. A quasi-equilibrium existsdue to the slowness of the energy relaxation via electron-ion collisions. Then a quasi-thermodynamic situation ex-ists for time scales τ such that τ ee ≪ τ ≪ τ ii ≪ τ ei ,where τ ss ′ is the temperature relaxation time [21] be-tween species s and s ′ . The NPA approach can be gener-alized to two-temperature WDMs and explicit NPA cal-culations are available [21, 22].Although the potentials in NPA are constructed in lin-ear response theory, they include the non-linear effectsbrought in via the Kohn-Sham calculation. Hence themethod is applicable to a wide range of densities andtemperatures. Here we examine simple fluids, mixturesof fluids, and complex fluids like liquid C, Si, containingtransient covalent bonds, from low T to high T . It isshown that these NPA pair-potentials work better thancommon multi-center potentials, and inexpensively re-cover pair-distribution functions and other properties ingood agreement with the best available DFT simulations,for well studies systems like Al, Be, Li, etc., and for com-plex fluids like C, Si and Ge. NPA calculations recoverthe properties of supercooled liquid Si at 2.57 g/cm justbelow the melting point of silicon, with only a fraction ofthe computation effort needed using VASP-type N -centercalculations [20].These methods, based on the two-component DFT ofions and electrons enable complex WDM calculations tobe reduced to simple one center Kohn-Sham calculationsthat require nothing more than a laptop. These reachthe accuracy of standard DFT calculations implementedby large codes like ABINIT [16] or VASP [17], whereadvanced gradient-corrected meta functionals have to beused because the N -center electron density is very com-plex and highly non-uniform.In the following we present our theory of pair poten-tials as follows. The NPA is introduced as a rigorousDFT concept and then various simplifications are indi-cated. We present a computational procedure for deter-mining the equilibrium one-body densities n ( r ) , ρ ( r ), andthe needed XC-functionals. Then we derive simple linear-response pseudo-potentials from the n ( r ) so obtained.These pseudopotentials apply to interactions mediatedby the continuum electrons, also called valance elec-trons, metallic electrons, or ‘free electrons’. The atomiccores are described by the bound states obtained fromthe NPA calculation. The pseudopotentials and bound-densities are used to construct bound-bound, bound-free,and free-free contributions to the pair potentials. Van derWaals type interactions occur in the bound-bound inter-actions. Then we compare our potentials with Stilllinger-Weber type multi-center potentials, and also with effec-tive medium theories. It is noted that, unlike the short-ranged potentials of multi-center empirical models, thelong-range pair-potentials of the NPA, with many min-ima that register with the second, third and higher neigh-bour effects in the pair-distribution functions, bring inmulti-ion correlations. An analysis of the potential feltby an ion in a fluid is used to show how these effects areincluded via the Onstein-Zernike equation and the mod-ified hyper-netted chain equation. Examples are given,including a comparison with the 40 parameter ‘glue po-tential’ of aluminum. II. DETAILS OF THE NPA MODEL
A brief description of the NPA model is given for con-venience and for defining the notation used. Hartreeatomic units ( | e | = ¯ h = m e = 1) will be used unlessotherwise specified. The systems explicitly studied arewarm-dense fluids and hence have spherical symmetry,although the NPA method applies equally well to low-temperature crystalline solids [23–25]. An application to surface ablation and hot phonons has been reported,e.g., [26] recently. However, the early NPA model (e.g., ofZiman) sought mainly to create a computationally conve-nient short-ranged object named the NPA, so that incon-venient long-ranged Coulomb fields of the ionic chargesare neutralized.The objective of the NPA model used here [2, 3, 25]is to implement a simplification of the full DFT modelfor electrons and ions given in Ref [27] and illustratedby an application to hydrogen. Its objectives have tobe distinguished from early NPA models, or from var-ious currently-available ion-sphere (IS) or average-atom(AA) models [28–31] used in many laboratories. Manyof them invoke simplifications which are outside DFT. Infact, the early formulations used an expansion of the ion-ion interactions in one-body, two-body, etc., terms andinvoked assumptions of non-overlapping electron distri-butions, superposition approximations, neglect of higherorder ionic and electronic correlations.They also invokedthe Born-Oppenheimer approximation. The formal the-ory does not need any of these assumptions [27]. How-ever, many of these approximations are invoked [2, 25]for numerical simplicity.Our NPA model is based on the variational property ofthe grand potential Ω([ n ] , [ ρ ]) as a functional of the one-body densities n ( r ) for electrons, and ρ ( r ) for ions. Onlya single nucleus of the material to be studied is used andtaken as the center of the coordinate system. The otherions (“field ions”) are replaced by their one-body densitydistribution ρ ( r ) as DFT asserts that the physics is solelygiven by the one-body distribution; i.e., we do not dis-card two-body, three-body, and such information as theyget included via exchange-correlation (XC)-functionals.Furthermore, in the case of ions, the exchange effects areusually negligible although the ion-ion manybody func-tional will also be called an XC-functional for generality.Hence unlike in N -center DFT codes like the VASP orABINIT, the NPA is a type of one-center DFT. Conse-quently, there is no highly inhomogeneous N -ion poten-tial energy surface, or a complicated electron distributionresident on such a multi-center potential energy surface.Hence a local-density approximation (LDA) to the elec-tron XC-functional can be sufficient. In our calculationswe have used the finite- T electron XC-functional by Per-rot and Dharma-wardana [32] in LDA, while most large-scale codes implement the T = 0 XC-functional inclusiveof generalized gradient corrections. The finite- T func-tional used is in good agreement with quantum Monte-Carlo XC-data [3] in the density and temperature regimesof interest.The artifice of using a nucleus at the origin convertsthe one body ion density ρ ( r ) and the electron density n ( r ) into effective two body densities in the sense that ρ ( r ) = ¯ ρg ii ( r ) , n ( r ) = ¯ ng ei ( r ) (1)The origin need not be at rest. However, most ions areheavy enough that the Born-Oppenheimer approxima-tion is valid. Here ¯ ρ, ¯ n are the average ion density andthe average free electron density respectively, and prevailfar away from the central nucleus. Any bound electronsare assumed to be firmly associated with each ionic nu-cleus and contained in their “ion cores” of radius r c suchthat r c < r ws , r ws = [3 / (4 π ¯ ρ )] / . (2)Here r ws is the Wigner-Seitz radius of the ions. The cor-responding electron sphere radius, based on ¯ n is denotedby r s . In some cases, e.g., for some transition metals,and for continuum resonances etc., this condition for acompact core may not be met, and additional steps areneeded. Here we assume that r c < r ws unless statedotherwise. The DFT variational equations used here are: δ Ω[ n, ρ ] δn = 0 (3) δ Ω[ n, ρ ] δρ = 0 . (4)These directly lead to two coupled KS equations wherethe unknown quantities are the XC-functional for theelectrons, the ion XC-functional, and the electron-ionXC functional [1]. If the Born-Oppenheimer approxima-tion is imposed, the electron-ion XC-functional may alsobe neglected. Approximations arise in modeling thesefunctionals and in simplifying the decoupling of the twoKS equations [2, 36]. The first equation gives the usualKohn-Sham equation for electrons moving in the exter-nal potential of the ions. This is the only DFT equationused in the N -center DFT-MD method implemented instandard codes like the VASP, where the ions define aperiodic solid whose structure is varied via N -atom clas-sical MD, followed by a Kohn-Sham solution at each step.In the NPA the one-body ion density ρ ( r ) rather than an N -center ion density is used. It was shown in [27] thatthe DFT equation for the one-body ion density can beidentified as a Boltzmann-like distribution of field ionsaround the central ion, distributed according to the ‘po-tential of mean force’ well known in the theory of fluids.Hence the ion-ion correlation functional F xcii was identi-fied to be exactly given by the sum of hyper-netted-chain(HNC) diagrams plus the bridge diagrams. The bridgediagrams cannot be evaluated exactly but may be mod-eled from a hard-sphere fluid [37], or evaluated using MD.The mean electron density ¯ n can also be specified asthe number of free electrons per ion, viz., the mean ion-ization state Z . It is also denoted as ¯ Z since a plasmacontaining a mixture of ionization states Z i + with com-position fractions x i will appear as a single average ionof charge ¯ Z = P i x i Z i + . This average description canbe sharpened by constructing individual NPA models foreach integer ionization and using them in the thermody-namic average [2]. Such a theory can yield the compo-sition fractions x s in the plasma, although the ABINITand VASP calculations cannot yield such information di-rectly.Although the material density ¯ ρ is specified, the meanfree electron density ¯ n is initially unknown at any given temperature, as it depends on the ionization balancewhich is controlled by the free energy minimization con-tained in the DFT Eq. 3. Hence a trial value for ¯ n (i.e,equivalently, a trial value for ¯ Z ) is assumed and the ther-modynamically consistent ρ ( r ) is determined. This isrepeated till the target mean ion density ¯ ρ is obtained.This also implies that NPA implements a rigorous calcu-lation of the ionization balance in the system. In effect,a full many-body extension of the Saha equation via itsDFT free-energy minimum property is implemented inthe NPA model to give the mean ionization needed inthe theory.The mean number of electrons per ion, viz., ¯ Z , is oftenreplaced by Z in this paper when no confusion arises.It is an experimentally measurable quantity. It is rou-tinely measured using Langmuir probes [38] in atmo-spheric and gas discharge plasmas, or from optical andother measurements of appropriate properties, e.g., theconductivity and the XRTS profile [39] for general WDMsystems. The argument that ‘there is no quantum me-chanical operator’ corresponding to ¯ Z , and hence it isnot an observable, is incorrect. The system studied isnot a pure quantum system but a system attached to aclassical heat bath. Here T, µ, ¯ Z are Lagrange multipli-ers associated with the conservation of energy, particlenumber, and charge conservation respectively. They donot correspond to simple quantum operators whose meanvalues are T, µ or ¯ Z , although complicated constructionscan be given. A. computational simplifications
The Kohn-Sham equation has to be solved for a singleelectron moving in the field of the central ion; its iondistribution ρ ( r ) = ¯ ρg ( r ) is modified at each iterationwith a corresponding modification of the trial ¯ Z untilthe target material density is found. However, it wasnoticed very early [2, 25] that the Kohn-Sham solutionwas quite insensitive to the details of the g ( r ). Hence asimplification was possible [2, 25, 27]. The simplificationwas to replace the trial g ( r ) at the trial ¯ Z by a cavity-likedistribution: g cav ( r ) = 0 , r ≤ r ws , g cav ( r ) = 1 , r > r ws (5) ρ cav ( r ) = ¯ ρg cav ( r ) , (6) V cav ( k ) = V k ¯ Z (3 /K ) [sin( K ) − K cos( K )] (7) V k = 4 π/k , K = kr ws . (8)Here the r ws is the trial value of the ion Wigner-Seitzradius, based on the trial ¯ n . Hence, adjusting the g cav ( r )at each iteration requires only adjusting the trial r ws toachieve self-consistency. The self-consistency in the iondistribution is rigorously controlled by the Friedel sum-rule for the phase shifts of the KS-electrons [27]. Thisensures that ¯ ρ = ¯ n/Z . Thus, a valuable result of thecalculation is the self-consistent value of the mean ion-ization state Z , which is both an atomic quantity and athermodynamic quantity. Here too we note that a sophis-ticated quantum many-body version of a Saha equationis also solved within this DFT calculation in determining¯ Z . The pseudopotentials, pair-potentials, the EOS etc.,depend directly on the accuracy ¯ Z .In early discussions of the neutral pseudoatom as ap-plied to solids [23, 24], no attempt was made to introducea DFT of the ions, and an analytically convenient poten-tial neutralizing the long-range Coulomb interaction ofthe central ion was added to produce a computationallyconvenient object with a short-ranged potential. In suchstudies, and in MD simulations of Coulomb potentials,a neutralizing Gaussian charge distribution is sometimesused. The spherical cavity potential used here is a betterapproximation to a typical ion-ion g ( r ).Here we note several simplification used in imple-menting the NPA. Given that the electron distribution n ( r ) obtained self consistently can be written as a core-electron (bound electron) term and a free-electron termwhen the core electrons are compactly inside the ionWigner-Seitz sphere, we have: n ( r ) = n c ( r ) + n f ( r ) , ∆ n f ( r ) = n f ( r ) − ¯ n (9) a f ( r ) = ∆ n f ( r ) (10)The core-electron density (made up of “bound electrons”)is denoted by n c ( r ). The free electron density n f ( r ) isthe response of an electron fluid (see sec. III C) to thecentral ion and to an ion-density deficit (a cavity) nearthe central ion that mimics ρ ( r ). It contributes to thepotential on the electrons. The response of a uniformelectron gas (i.e., without the cavity) to the central ioncan be obtained by subtracting out the effect of the cavityusing the known static interacting linear response func-tion χ ( k, ¯ n, T ) of the electron fluid. That is, from nowon we take it that the charge density n f ( r ) and the den-sity pile up ∆ n f ( r ) are both corrected for the presenceof the cavity. For simplicity we use the same symbols asbefore when no confusion arises. Furthermore, in dealingwith atoms ‘a’, ‘b’ etc, we refer to the ‘density pileup’or displaced free electron density ∆ n f ( r ) as a f ( r ), withits Fourier transform denoted by a f ( k ). In Fig. 1(a) wedisplay the calculated n f ( r ) for an Aluminum ion in ahot plasma.A key difference between many average-atom modelsand the NPA is that the free electrons are not confinedto the Wigner-Seitz sphere, but move in all of space rep-resented by a very large ‘correlation sphere’ of radius R c .This is ten to twenty times the Wigner-Seitz radius of thecentral ion [2, 29]. This ensures that the small- k limit ofthe structure factor S ( k ) is accurately obtained. III. PSEUDOPOTENTIALS ANDPAIR-POTENTIALS FROM THE NPA MODEL.
As the NPA reduces the multi-center WDM to a sin-gle ion and its environment described by n ( r ) , ρ ( r ), pairpotentials must be constructed from the appropriate one r [a.u.] n f (r) [ a . u . ] free electron pileup k/k F -1-0.500.5 F o r m f ac t o r M ( k ) M(k)= ∆ n f (k)/|ZV k χ( k)| n Al, 2.7 g/cm T=10 eV r ws (a) (b) FIG. 1. (Online color) Panel (a). The free electron density n f ( r ) around an Al nucleus placed in its electron fluid ofaverage density ¯ n corresponding to ¯ Z = 3, ρ = 2 . , at10 eV. The bound electron core is well inside r ws , and n f ( r )develops oscillations as the free-electron states are orthogonalto core states. Panel (b). The k - dependent pseudopotentialform factor M ( k ), calculated from the electron pileup ∆ n f ( k )and the electron response function χ ( k, T ). Here V k is theCoulomb potential 4 π/k . The beahviour for large k/k F > k F is the Fermi wave vector. body distributions, viz, the core-electron density n c ( r ),and the ‘free’ electron density n f ( r ).The pseudopotential U ei ( r ), or its Fourier transformedform U ei ( k ) enables a rigorous and useful separation ofthe contributions of bound and free charge densities tothe pair potential which is a two-center property. Weconstructed it from one-center results to avoid two-centercomputations. Thus the pseudopotential is constructedto be a linear response property where possible. In fact,unlike in crystals, the spherical symmetry of fluids allowsone to use simple local pseudopotentials, without havingto deal with angular-momentum dependent (i.e., nonlo-cal) pseudopotentials. These potentials are analogous tothe s -wave potentials used in the N -center DFT codes,but differ from them in that NPA pseudopotentials areconstructed to be weak potentials, allowing the use ofsecond-order perturbation theory in their manipulation.Furthermore, unlike non-linear pseudopotentials usedin ABINIT or VASP, these pseudopotentials are ρ, T de-pendent, and are constructed in situ during the calcula-tion. The range of validity of such linearized pseudopo-tentials has been discussed elsewhere [40].A rigorous discussion of pair-potentials U ab betweentwo atoms of species ‘a’, and ‘b’, including the effect ofcore electrons is given in the Appendix B of Ref. [2].Core-core interactions are important for atoms includingand beyond argon, viz., Na, K, Au, etc., with large coresand zero or low ¯ Z . In contrast, core-core interactions inhigh Z ions like Al , Si at normal density are small incomparison to the interactions via the metallic electrons.In dealing with core-core interactions we illustrate theconcepts used, by a discussion of the case of the weaklyionized argon WDM system. We denote the core electrondistribution n c ( r ) of the atom ‘a’ by a c , while a f denotesthe displaced free electron distribution, δn f ( r ), for theatomic center of type ‘a’ (and similarly for ‘b’). Thetotal pair-interaction V ab is of the form: V ab = V ( a c , b c ) + V ( a c , b f ) + V ( a f , b c ) + V ( a f , b f )= V cc ( r ) + V cf ( r ) + V ff ( r ) (11)We discuss each of these terms separately. A. The interaction mediated by free electrons
The last term of Eq. 11, viz., V ff ( r ) is the familiarion-ion interaction mediated by metallic electrons. Thisis adequately evaluated in second-order perturbation the-ory when the interactions are weakened due to free-electron screening effects. The linearized electron-ionpseudopotential described below can be used for WDMsystems and liquid metals unless the free electron den-sity and the temperature ( T /E F ) are low, when linear-response methods become invalid. In fact, when the pair-potential develops a negative region significantly largerthan the thermal energy, permanent chemical bonds areformed; then the liner-response methods used here can-not be used.The electron-ion pseudopotential U ei a ( k ) of the ion ‘a’ isevaluated in k -space via the displaced free-electron den-sity a f ( k ), and the electron response function χ ( k, T ).The pseudopotential is: U ei a ( k ) = a f ( k ) /χ ( k, T ) (12)= − Z a V k M ( k ) , V k = 4 π/k (13)Here M ( k ) is the form factor of the pseudopotential (seeFig. 1). The NPA calculation via Eq. 12 automaticallyyields the pseudopotential inclusive of a form factor. Thismay be fitted to a Heine-Abarankov form which fits a coreradius r c and a constant core potential V c . Alternatively,a more complex parameterized form may be needed, asin Ref. [40]. However, using the numerical tabulation of U ei a ( k ) is more accurate and avoids the fitting step.The fully interacting linear response function χ ( k, T )of the uniform electron fluid at the temperature T is dis-cussed in sec. III C. So, for identical ions ‘a’, we have theform: V ff aa ( k ) = Z a V k + | U ei a ( k ) | χ ( k, T ) (14)= Z a V k + | a f ( k ) | /χ ( k, T ) . (15)In the case of two different atomic species, we readilyhave: V ff ab ( k ) = Z a Z b V k + (16) | U a ( k ) U b ( k ) | χ ( k, T ) (17)In the NPA theory for mixtures, Z a and Z b are inte-gers, while in the simple average-atom form of the NPA, r/r s V ii (r) / T & g (r) V(r) T=0.5 eVg(r) T=0.5 eVg(r) T=2.0 eVV(r) T=2.0 eV
H density ρ = 1.0 g/cm r s = 1.392 FIG. 2. (Online color)The proton-proton pair potentials andPDFs at 0.5 eV and 2 eV in a Hydrogen plasma at r s =1.391,with ρ =1.0 g/cm , at T significantly higher than a possiblephase transition and contains protons and electrons.The pair-potentials contain Friedel oscillations whose minima line upwith the peaks in the g ( r ). the ¯ Z s = Σ s x s Z s , is used to represent the mixture withcomposition fractions x s via a single calculation. In amixture calculation, separate NPA calculations for eachcomponent are needed.The simplest pseudopotential, valid for point ions, r c = 0, uses a form factor of unity. The long-wavelengthform ( k →
0) of the response function applies far away( r ≫ r ws ) from an ion in the uniform-density region of aplasma. Then it depends only on the screening wavevec-tor k s ( T ). Given such approximations, the electron re-sponse function χ ( k, T ) and the pair-potential reduce tovery simple forms. V k χ ( k ) = k s /k , V ff ab = 4 πZ a Z b / ( k + k s ) (18) k s = k se ( T ) = 4 πT Z ∞ k dkn ( k ) { − n ( k ) } (19) n ( k ) = 1 / (cid:2) { ( k / − µ ( T )) /T } (cid:3) (20) V ff ab ( r ) = ( Z a Z b /r ) exp( − k s r ) (21)Here k s = k se is an electron screening wave vector in thefluid at the temperature T and at the electron chemicalpotential µ ( T ). This reduces to the Thomas-Fermi valueat T →
0, and to the Debye-H¨ukel value as T → ∞ .However, most scattering processes for T /E F ≤ k F within a thermal win-dow of the Fermi energy. Hence the use of the k → T .The real space form of the long-wavelength screen-ing formula is well known as the “Yukawa potential”, Z a Z b exp( − k s r ) /r . It has been used in particle physics,physical chemistry and in plasma applications [42–44].The point ion model is unsatisfactory even for protonsfor T < E F as discussed below. The proton-proton pair-potential is a most demanding case for the NPA methodbecause the charge pile up around the proton is highlynon-linear, and the methods used here become inappli-cable without further modification when very few freeelectrons are present, while partial ionizations strictlyrequire the use of a mixture of ions of different ioniza-tions. Forms like H − , H, H + as well as quasi-moleculartransient forms like H +2 , may occur in H plasmas neara phase transition. The fully ionized case (i.e., a pro-ton) has no bound core, and yet the point-ion modelfails since a form factor associated with the non-linearcharge pile up is needed. Nevertheless, we find that themethod works with adequate accuracy, and even picks upthe effects of transient H +2 in the medium when relevant.we give an illustrative example below, using fully ionizedhydrogen.Fig. 2 displays the the H-H pair potential and PDFsat 0.5 eV and 2 eV for a hydrogen plasma at r s ≃ . ρ = 1 . . The plasma is fully ionized, althoughclose to a first-order transition from a molecular liquidto a conducting atomic fluid, believed to occur near T ∼ g ( r ) in Fig. 2 from 1.5 < r/r s < . +2 states as well as withH + packing effects. The peak due to purely packing ef-fects occur near r/r ws ∼ .
6. The nominal H +2 bondlength of ≃ T = 0 is extended ina plasma due to screening and temperature effects. Thebond becomes transient, with a broader range of lengths.It should be noted that the NPA pair-potentials arevery long-ranged, encompassing some tens of Wigner-Seitz radii. The aluminum pair-potential near its melt-ing point requires using a pair-potential extending overat least a dozen atomic shells if the phonon spectra are tobe accurately recovered [26]. This should be compared tothose used in the pair-potential part of empirical mod-els like the SW potential which extends to only aboutone Wigner-Seitz radius. Furthermore, the minima inthe Friedel oscillations of the NPA potentials correspondto shells of atoms where the third, fourth and higherneighbours are positioned, as seen by a comparison withthe PDFs generated from these potentials. The PDFsare generated using the NPA potentials in classical MD,or in a hyper-netted-chain equation inclusive of a bridgeterm.A potential with wider applicability to uniform flu-ids than the Yukawa model is obtained using the ran-dom phase approximation (RPA) to the response func-tion. This retains the all important Friedel oscillationsin the pair potential for T /E F <
1. At high temper-atures this reduces to the Vlasov approximation. The
D-D distance, r/r s g dd (r) g(r,T=0.75 eV)+0.75g(r,T=0.5 eV)+0.50g(r, T=1.0 eV)+0.25g(r,T=2.0 eV) D-D Molc. pre-peak
Deuterium, r s =2 a.u. ρ= FIG. 3. (Online color)The deuteron-deuteron pair distribu-tion functions display the onset of transient D-D molecularforms (e.g., D +2 ) as the temperature is lowered from 2 eV to0.3 eV, for a deuterium plasma at r s = 2. The D + are fullyionized and only transient molecular forms exist, while thereare no atomic bound states. pair-potential, the detailed computations of the pair dis-tribution functions (PDFs), free energies and EOS prop-erties of a plasma of point charges interacting via theRPA screened Coulomb potential is an important refer-ence point. These have been computed as a functionof the plasma parameter Γ = Z/ ( r ws T ) by Perrot [45].However, the assumption of negligible core size (point-ion approximation) is not adequate for most ions exceptat high T . Neither the RPA form, nor the Yukawa formsatisfies the compressibility sumrule.A complete treatment of ion-ion interactions needscore-core interactions as well as the inclusion of a formfactor in the electron-ion pseudopotential U ei a ( k ). We useEq. 12, which is an accurate yet simple pseudopotentialdirectly available from any average-atom calculation. Wehave successfully used the NPA U ei ( k ) in many diversematerials and WDM systems ranging from solid Be, Na,Al, to plasma states of Be, Li, C, Ga, Al, or Si even as asupercooled liquid.The pseudopotential resulting from Eq. 12 is a simplelocal potential ( s -wave potential). If the density displace-ment a f ( k ) is expanded in terms of the contributionsfrom different angular momentum states l of the Kohn-Sham eigenstates, then a non-local ( l -dependent) pseu-dopotential can be constructed. Such non-local forms areneeded in solid-state applications where spherical sym-metry is lacking. However, the simple s -wave form givenin Eq. 12 seems to work well for fluids and their EOS andstatic transport coefficients.Aluminum is regarded as a ’difficult’ material by thosewho develop effective medium potentials for it. If Fig. 4we display Al-Al pair potentials for several tempera-tures and compressions, including for a case where theions are held at 1 eV, while the electrons are at 1 eV, r [Å] V ii (r)[ e V ] T e = 1 eVT e = 5 eVT e = 10 eV r/r ws V ii (r) / T & g (r) V( ii r) g(r) Al, 2.7 g/cm (a) (b) T i = 1 eV T=T i =T e =1 eV T=0 XC
Al, 2.7 g/cm FIG. 4. (Online color) (a)The Al-Al pair potentials at nor-mal density from the NPA for WDMs in thermal equilibriumat T = 1 eV, and two WDM states with T e = T i . The Fermienergy E F ∼
12 eV, and the electrons in this system are par-tially degenerate. If the electron XC is treating using the T = 0, as is done in VASP and ABINIT calculations, thereis a slight error in the depth of the first energy minimum.(b) The Al pair-potential at 1 eV (black line), as well as thecorresponding g(r) shown as a red line with triangles. Theenergy is scaled by T , and the position by r ws , unlike in theleft panel. The first maximum in the g(r) is at a positive en-ergy and pulled inwards from the first minimum in the pairpotential, and yet it is the preferred first-neighbour position.This is because the denser fluid is favoured by the XC-energyof the electron fluid. The second and third neighbours, corre-sponding to the secondary peaks in the PDF are more closerto the minima of the Friedel oscillations. or at 5 eV and 10 eV. The potentials show how theFriedel oscillations disappear at high electron temper-atures. These minima, which determine the locationof nearest-neighbour, next-nearest-neighbour positionsbring in the multi-center features that empirical modelslike the SW potential put in by hand. However, Fig. 4shows that the peaks and troughs in the g ( r ) not entirelydetermined by the pair potential as the ‘volume energy’of the electron fluid plays a part, and largely determinesthe position of the first peak in a system like liquid alu-minum or liquid carbon (see discussion in Ref. [34]) wherethe electron density is high.The PDFs g ( r ) can be obtained from the pair-potentialeither using classical MD, or using the modified HNCequation (MHNC). The g ( r ) shown in Fig. 4(b) has beenobtained using the MHNC and a Bridge term based ona hard sphere liquid with a packing fraction of 0.2996.It is also of interest to see how well these NPA PDFsof one-center DFT agree with direct simulations usingan N -center DFT procedure like that of the VASP orABINIT. In Fig. 5 we give comparisons of NPA for twoPDFs with those of Recoules et al. [46], with N =64 atomsin the simulation cell. The PDFs obtained from DFT-MD simulations are sensitive to the type of electron XC- r [a.u.] g (r) DFT-MD, 0.86 eV (10,000 K)NPA-MHNC, 0.86 eV DFT-MD, 2.58 eV (30,000 K)NPA-HNC, 2.58 eV
DFT-MD by Recoules et al.Al, 2.7 g/cm FIG. 5. (Online color)The NPA g ( r ) for aluminum at twotemperatures are compared with the DFT-MD simulationsof Recoules et al [46]. The simulations use 64 atoms andthe PBE functional. Better agreement near the first peak isobtained from simulations with more atoms. functional used, and a higher N is needed to get betterstatistics. However, the NPA g(r) and the DFT-MD canbe considered to be in good agreement. B. Contributions to the potential from thecore-electron density
Equation 11 contains V cc and V cf where the core-electron density contributes to the pair-potential. In thefollowing we show that the neglect of core-core interac-tions can lead to very serious differences, e.g., speciallyin weakly ionized plasmas containing neutral species.The term V cc ( r ) contains V ( a c , b c ) which is the in-teraction between the density distributions of the boundelectrons in the atom ‘a’ with those in ‘b’. These are notthe unperturbed NPA densities n c ( r ) defined in Eq. 9,but the densities that result from the perturbation ofeach density by the other, in the plasma environment.The unperturbed core density of the ion ‘a’ is denotedby n c ( ~r a ) , ~r a = ~r − ~R a and similarly for the ion ‘b’, with R = | ~R a − ~R b | defining the scalar separation of the twoatoms. The contribution to the pair-potential is: V cc ( R ) = Z a Z b /R + Z dr dr ′ n c ( ~r a ) n c ( ~r ′ b ) | ~r − ~r ′ | +∆ ES ( R ) + ∆ xc ( R ) (22)The potential contains an electrostatic correction∆ ES ( R ) as the core densities are modified by the inter-action. Instead of doing a two-center calculation for thisterm, second-order perturbation theory based on the po-larisability of the core density is sufficient specially whenthere are free-electron screening effects, as in plasmasor metals. An electron XC-term ∆ xc ( R ) is associatedwith the density modification. This includes a correctionto the kinetic energy functional as well as terms arisingfrom XC-effects. That is, evaluating ∆ xc ( R ) accuratelyis somewhat complicated. A strictly density-functionalapproach is discussed in Appendix B of ref. [2].A simplified approximate procedure is as follows.The core-core interaction can be written as a sum ofmonopole, dipole, quadrupole and higher terms. In sys-tems with free electrons (e.g., plasmas), only the dipoleterm of the expansion is of importance due to screen-ing effects. The frequency-dependent part of the dipoleinteractions brings in the van der Waals (vW) type ofcontributions (“dispersion forces”). These have been dis-cussed extensively within DFT [47, 48]. While they areeasily included in the NPA approach, they are hard toinclude in the usual N -center DFT used in codes like theVASP because of their strong non-locality. One methodused in N -center DFT calculations is to decompose the N -center charge density into N individual charge distri-butions (each like an NPA). Using maximally localizedWannier functions on each atomic center and a polariz-ability analysis for each center is one typical approach.The need for such procedures is an artifact of the usual N -center DFT approach where even a fluid is representedas an average over an ensemble of MD generated solids.In the NPA scheme we already have such single-centerdistributions in the appropriate form and dipole forcesand vW effects are easily included when needed.In practical WDM calculations where free electrons arepresent the metallic interaction V ff ab ( k ), Eq. 16, is domi-nant. Then a fairly simple procedure may be sufficient forincluding core-core interactions when they become rele-vant. For instance, given an atom ‘a’ with an argon-likecore, and an atom ‘b’ with a krypton-like core, we use theknown Ar-Kr interaction given as a multipole expansioninclusive of dispersion forces, and screen them using thefree-electron dielectric function of the WDM electrons atthe appropriate free-electron density ¯ n and T .We illustrated this using an example, namely, a mix-ture of neutral argon atoms Ar , together with argon ionsAr Z i + with appropriate integer ionizations Z i [49]. In anargon plasma with, say, ¯ Z = 0 . ∼ ∼
30% are singly ionized Ar + . Treating such asystem needs the neutral Ar-Ar interactions which arepurely core-core interactions. We also have Ar-Ar + aswell as Ar + -Ar + interactions, i.e., three pair-potentialsand three ion-ion PDFs associated with them. A singleaverage-atom NPA is inadequate for accurately estimat-ing physical properties of such a system.In evaluating physical properties of the Ar, Ar + mix-ture, say, the self-diffusion coefficient, there are two self-diffusion coefficients, and one inter-diffusion coefficientas well. They are also constrained by the compressibilitysum rule. In such instances, the meaning of the self-diffusion coefficient obtained from an N -center DFT cal-culation (e.g., a VASP calculation) needs to be examinedmore carefully as no information about ionization states r [a.u.] V (r) / T neutral Ar-ArAr-Ar + Ar + -Ar + r [a.u.] V (r) / T neutral Ar-ArZ=0.3, no-core interactionsT = 2 eV, 1.395g/cm FIG. 6. (Online color)The screened pair-potentials Ar-Ar,Ar-Ar + Ar + -Ar + of a two-component mixture of argon at 2eV. The single NPA pair-potential V ff a , a ( r ), a=Ar . is ob-tained (thick magenta curve) for ions with ¯ Z =0.3, and doesnot include core-core interactions of individual ions is available in such calculations.Evaluations of core-core potentials for neutral argonatoms in a medium without free electrons give resultsclose to parameterized potentials similar to the Lennard-Jones(LJ) or more advanced rare-gas potentials. Thepresence of free electrons leads to screening of the core-core interaction. Thus, at the LJ-level of approxima-tion, V cc ( k ) for two neutrals immersed in the appropri-ate electron gas is approximately given by V LJ ( k ) { V k χ ( k, T ) } .For two charged ions, the major interaction is via freeelectrons, and is given by the usual form, Eq. 16. The ma-jor part of the interaction between the ion and the neutralatom is the energy of the induced dipole of the neutralatom interacting with the electric field of the ion, mod-erated by screening. There is also a shell-shell repulsion,giving rise to a modified LJ-like potential. The polariz-ability and other parameters can be evaluated using thebound electron densities of the Ar atom and the Ar + ionobtained from the respective NPA calculations, or usingknown LJ parametrizations. Corresponding calculationsusing the VASP code etc., and using NPA results, havebeen reported by Stanek et al., Ref. [49]. In Fig. 6 wedisplay the three pair-potentials relevant to argon at nor-mal density and T = 2 eV, when ¯ Z = 0 .
3, and treated asa mixture. However, if the system is treated as a singlecomponent fluid of average atoms with a charge of ¯ Z =0.3,the resulting pair-potential evaluated without includingcore-core corrections is found to be quite different (rightpanel, Fig. 6) from those inclusive of core corrections.0 C. The response function of the uniform electronfluid
The response function is itself a functional of theone-body electron density and can be evaluated self-consistently within the calculation by solving the NPA-Kohn Sham equation for just the central cavity but with-out a nucleus. Then the free electron density pile up a f ( k ) is entirely due to the cavity potential V cav ( k ) whichis the electrostatic potential of a weak spherical cavity.Then Eq. 12 can be used in an inverse sense to calcu-late χ ( q, T ) in situ . While such an approach is useful asa control calculation, the following direct procedure hasbeen implemented in our codes.The interacting electron gas response function used inthese calculations includes a local-field factor chosen tosatisfy the finite temperature electron-gas compressibilitysum rule. χ ( k, T ) = χ ee ( k, T e ) = χ ( k, T e )1 − V k (1 − G k ) χ ( k, T e ) , (23) G k = (1 − κ /κ )( k/k TF ); V k = 4 π/k , (24) k TF = { / ( παr s ) } / ; α = (4 / π ) / , (25) V ii ( k ) = Z V k + | U ei ( k ) | χ ee ( k, T e ) . (26)Here χ is the finite- T Lindhard function, V k is the bareCoulomb potential, and G k is a local-field correction(LFC). The finite- T compressibility sum rule for elec-trons is satisfied since κ and κ are the non-interactingand interacting electron compressibilities respectively, atthe temperature T , with κ matched to the F xc ( T ) used inthe KS calculation. In Eq. 25, k TF appearing in the LFCis the Thomas-Fermi wavevector. We use a G k evaluatedat k → k instead of the more general k -dependentform (e.g., Eq. 50 in Ref. [32]) since the k -dispersion in G k has negligible effect for the WDMs of this study. IV. THE NPA PAIR-POTENTIAL APPROACHCOMPARED TO METHODS BASED ONMULTI-CENTER POTENTIALS
The pair-potential V ab ( R ) is normally understood asthe energy of the system with two ‘atoms’, or ‘ions’ asthe case may be (but referred to here as ‘atoms’) anddenoted by ‘a’ and ‘b’ and held at a separation R , com-pared to the limit R → ∞ when there is no interactionbetween them. This has a clear meaning if the system isin a vacuum and at sufficiently low T such that ‘a’ and‘b’ remain as compact objects. Thus the ‘bond energy’ oftwo hydrogen atoms calculated using standard quantumchemistry methods, at T = 0, has no ambiguity. How-ever, in an atomic or molecular fluid, or in a WDM, themedium contains other atoms and even free electrons,with T = 0.So one needs to include the effect of neighbouringatoms that are in the medium. This is exactly the prob-lem systematically addressed by many-body theories like DFT. The effect of the medium is a functional of the one-body densities of the components that make up themedium where the two atoms ‘a’ and ‘b’ are placed in.Instead of constructing the necessary functionals in asystematic way, keeping in mind that they should beone-body functionals, semi-empirical potentials (fittedto data bases etc.) have deployed multi-center modelsmimicking bonds, bond angles, their torsional propertiesand so forth, based on a preconceived ‘chemical bond-ing’ picture. If this ‘fitting’ had been directed to con-structing just the electron XC-functional, then we havethe effort found in quantum chemistry, colourfully de-scribed as constructing a “Jacob’s Ladder”, where in-creasingly more sophisticated electron XC-functionals areconstructed for electronic systems placed in an externalpotential of a finite number of nuclei, usually held fixed.In the following we illustrate the DFT-NPA approachof using pair-potentials and XC-functionals by comparingit with the Stillinger-Weber (SW) model as a well-knownexample of a multi-center potential useful for tetrahedralmaterials, be they in solid, liquid or even WDM stateslike molten silicon or liquid carbon. A. The N -body energy in the Stillinger Webermodel We briefly recount the details of the SW model for theconvenience of the reader. The quantity they model is thepotential-energy function φ for N identical interactingparticles. SW write the potential for Si as: φ (1 , , . . . , N ) sw = X i 959 a.u., a = 1 . 80 and ǫ = 2 . bond length inthe solid. Such a deep feature is not found in DFT-VASPor in the DFT-NPA potentials. Unlike in SW contain-ing only classical ions, the NPA and VASP calculationstreat silicon as a two-component system with an ion sub-system and an electron subsystem. The stability of thesystem arises from the ionic interactions as well as a vol-ume energy associated with the electron fluid. Hencethe DFT models have a comparatively shallow ion-ionpotential, while the SW potential has to be deeper toinclude the energy of the electron subsystem as a “bondenergy” component. Furthermore, while a very short-ranged pair-potential alone is insufficient to ‘stabilize’ theone-component tetrahedral structure, a long-ranged pair-potential inclusive of the ‘electron gas’ contribution tothe total energy covers all the possible electron-ion struc-tures that are quantum mechanically possible. As al-ready noted in discussing the potentials and PDFs of Al,Ar, hydrogen, and in previous publications on C, Si [34],Al [35], the energy minima in the Friedel oscillations ofthe ion-ion pair potential, together with the volume de-pendent energy effects of the electron fluid, were seen tocorrelate closely with the positions of the higher-orderneighbours beyond the nearest neighbour.A comparison of the Si-Si g ( r ) obtained from vari-ous calculations is given in Fig. 8. The PDFs fromNPA calculations, DFT-MD calculations using the SCANXC-functional, as well as the PBE [51] XC-functional r [Å] -2-101 V ii (r) [ e V ] -1200K, Z=4 NPA2.53 g/cm -1200K, S-W2.27 g/cm -1200K, DFT-MD1.50 g/cm -1200K, Z=3 NPA S i + bond Si FIG. 7. (color online) The pair-potentials from the NPA,from a VASP calculation with 216 atoms using the SCANXC-functional for supercooled liquid silicon at 1200K, and atthe metastable low density of 2.27g/cm are shown. The nor-mal melt density is higher, being at 2.57 g/cm prevails above1683K. The pair-potential for the unstable higher density liq-uid from NPA at 1200K is also shown. The Stillinger-Weberpair-potential for liquid Si at 1200K is contrasted with theother potentials. We also display the NPA potential for ex-panded molten Si at 1.5g/cm where a low-ionization stateSi occurs. The Si-Si nearest-neighbour “bond” distance isshortened for Si + ions, compared to the Si pairs. are shown, together with a typical classical MD simu-lation for the SW potential [52]. The one-center elec-tron density distribution used in the NPA is very sim-ple and smooth compared to the N -center n ( r ) usedin the VASP calculation. Hence the local density ap-proximation to exchange-correlation works efficiently andaccurately. The NPA calculation with hyper-netted-chain (HNC) equation is computationally much fasterand cheaper than the SW simulation, while also beinga first-principles DFT calculation. B. The ion-ion N -body corrections and the NPAmodel It is of interest to see how the pair-potential used inthe NPA brings in the three-body energy and such ‘multi-center’ terms via the ion-ion XC functional used in theNPA approach. Of course, if classical MD is used withthe NPA pair-potentials [22, 49, 53], then the ion-ion cor-relations are automatically built up by the long-rangepair-interactions. Such simulations yield the g ii ( r ), butthe calculation of the total energy requires contributionsfrom the electron subsystem in calculating the total en-ergy, as discussed in Ref. [2]. In using the SW potential,the g ii ( r ) and the SW-potential are sufficient to deter-mine the total energy as there is no electron subsystem.Here we consider the question of many-ion effects using2 r [Ang] g (r) NPA-LDA-HNCDFT-PBE-MDDFT-SCAN, Remsing et al.Stillinger-WeberLiquid Si. S i - S i bond . A . N - N - N , . A o o FIG. 8. (Online color) The PDFs for molten Si slightly abovethe melting point calculated by several methods is displayed.Two DFT-MD calculations using the VASP, code, and withthe SCAN and PBE XC-functionals show sensitivity to thetreatment of electron exchange correlation used in the 216atom calculation. The NPA, using only one Si atom has amuch smoother electron distribution and obtains comparableresults with an LDA XC-functional. The indicated nearest-neighbour (N-N) and N-N-N distances are for solid Si. TheS-W model shows a spurious hump near 3-3.5 ˚A . the theory of integral equations for ion correlations.In contrast to the SW-potential, the NPA pair poten-tials do not include any pre-assigned structural charac-teristics. Such pre-assigned structures are usually validonly in a specific regime of ionization (e.g., carbon witha valance of four), density and temperature. The pair-potentials and the associated XC-functionals of the NPAapproach generate the appropriate lowest energy struc-ture where the local ordering may be tetrahedral, face-centered etc., or a thermodynamic mixture of manystructures with no predominant structure.We consider a uniform fluid consisting of atoms of onetype for simplicity. The DFT Eq. 4 for the ion subsytemsolves for the minimum energy ion distribution ρ ( r ) andgives the form: ρ ( r ) = ¯ ρg ( r ) = ¯ ρ exp {− V I ( r ) /T } (35)where V I ( r ) is the density-functional potential felt by anion I at the radial location r , in the presence of an ionat the origin. This potential can be written within DFTas: V ( r ) = V I ( r ) + V (mean field) + V xcii (r) (36) V xcii = V cii ( r ) = δF xcii [ ρ ( r ) , n ( r )] /δρ ( r ) (37)The first term on the r.h.s. is the ion-ion pair poten-tial between the central ion and the ion I at r , viz.,the NPA pair-potential. Any arbitrary ion is denotedby the lower-case i . The ion I is also subject to the self-consistent average potential at r from all the field ions in the medium, i.e., from the density ρ ( r ) = ¯ ρg ( r ). Theeffect of the electron distribution n ( r ) has to be includedin the mean field potential. The electron density is givenby the simultaneous solution of the KS-Eq. 3. The meanfield potential is just the Poisson potential Z d~r ′ ( ¯ Zρ ( ~r ′ ) − n ( ~r ′ )) / | ~r − ~r ′ | when evaluated at the level of a monopole expansion, i.e.,neglecting core effects and pseudopotential form factors,for simplicity.Terms not contained in the mean-field are in the ionic-XC term V xcii ( r ). We also neglect electron-ion XC effects(see Ref. [1]). Since ions are classical, there is no ex-change, and this is purely a correlation term that can becalculated using classical statistical mechanics. The ioncorrelation free energy is formally given as a coupling-constant integration on the g ii ( r ) [2], just as for the elec-tron XC-functional. But g ii ( r ) is initially unknown. Useof an LDA V cii ( r ) based on the classical correlation en-ergy of the uniform classical Coulomb fluid, as done forelectron XC functionals turns out to be a failure, as re-ported in Ref. [27]. It is found that V cii ( r ) is highly non-local and even gradient expansions fail. Hence we exploitthe relationship of XC-functional with the correspondingpair-distribution functions as follows.To include the effect of the three-body and higher cor-relations neglected in the mean-field potential, we use theOrnstein-Zernike relation which expresses the total cor-relation function h ( r ) = g ( r ) − c ( r ). h ( ~r ij ) = c ( ~r ij ) + Z d~r k c ( ~r ik ) ρ ( ~r k ) h ( ~r kj ) (38)This equation is a simple algebraic equation in k -spacefor fluids of uniform-density. It brings in the interactionof the two atoms i, j with all possible ‘third atoms’ atthe location ~r k . The explicit correlation corrections areselected to be the sum of hyper-netted-chain (HNC) dia-grams excluding mean-field terms that have already beentaken into account in Eq. 36. That is, noting that r = | ~r ij | in the uniform case, − V cii ( r ) /T = ¯ ρ Z { h ( ~r ′ ) + V I ( r ) /T } h ( ~r − ~r ′ ) d~r ′ (39)Correlation terms not captured by the HNC sum are noteasily evaluated and are bundled into a ‘bridge diagram’term that is evaluated using hard-sphere models [37].In actual NPA calculations where a large correlationsphere of radius R c is used, all quantities are referred tothe limit r → R c . This is true for all XC-potentials aswell as other quantities like chemical potentials. That is, V xcss ′ ( r ) are effective values V xcss ′ − V xcss ′ ( R c ).3 C. Expression for total Free energy in the NPAmodel For comparison with expressions to be given below foreffective medium theories, we discuss the total Helmholtzfree energy F of a fluid of neutral pesudoatoms. The eval-uation of the free energy is more complicated than thatof the internal energy E tot , as calculations of the entropycontribution to F in DFT is not conceptually straightfor-ward especially in dealing with atomic bound states withpartial occupancies. The total free energy per atom for agiven ionic configuration (expressed as the density distri-bution ρ ( r ) = ¯ ρg ( r ) around a given nucleus) is schemat-ically summarized below, while the full expressions maybe found in Refs. [25] and [2]. F = F id + F + F em + F − E at (40) F = ¯ Zf h (¯ n, T ) = ¯ Z ( f k.e + f xc ) (41) F em = ∆ F + ∆ V ei (∆ n ) + ∆ V ee (∆ n ∆ n ′ ) (42) F = (1 / ′ X j V ( R j ) + cavity corrections (43)The energy of an isolated atom at T = 0 in its groundstate, E at , is used as the reference energy. Here F id isthe classical ideal-gas energy of the non-interacting ionsubsystem per atom. Each atom contributes ¯ Z free elec-trons to the electron fluid. F is the free energy of ¯ Z electrons in a uniform electron gas, with a kinetic energycomponent and an XC-component. The “embedding freeenergy” of a neutral pseudo atom is F em . This is definedas the difference between the free energy of the electronfluid of mean density ¯ n with and without the NPA. Itconsists of a correction to F , the Coulomb interactionof the pseudo-ion with the displaced electrons, as well asthe correction to the electron-electron repulsion energydue to the displaced electrons.The F term contains the free energy contributedthrough the ion-ion pair-distribution function and thepair-potential.The many-atom correlations dependent onthe ‘environment’ are included here via the ion-ion PDF.Hence this is effectively the pair-energy and the correla-tion corrections to the free energy of the classical ions.That is, in more conventional language, the ‘bonding en-ergy’ of two ions immersed in the fluid with average den-sities ¯ n and ¯ ρ , inclusive of the effect of their higher-orderneighbours brought in via the secondary peaks of thePDFs.While F is a smooth function of the electron density,sharp jumps may occur in F , and to a lesser extentin F em , signaling the onset of ionization changes, struc-tural changes or phase transitions. Example of such dis-continuities at liquid-liquid phase transitions (LPT) maybe found in Ref. [20] for liquid Si, while similar tran-sitions have been noted in theoretical studies of liquidcarbon as well [6, 7]. Interestingly, a molten-carbon LPTwas predicted in Ref. [6] using an empirical bond-ordermulti-center carbon potential, but this was not confirmed when DFT-MD simulations were done using the PBEfunctional [8] V. EFFECTIVE MEDIUM POTENTIALS ANDNPA POTENTIALS Unlike the purely classical SW model and effec-tive medium approach (EMA), embedded-atom models(EAM) pay attention to the existence of an electron sub-system, and begin by using approximations to T = 0DFT to include the “effect of the local environment” ofan ion in the medium. The Finnis-Sinclair potentials [55]used in metallurgical applications is based on a second-moment approximation to tight-binding but reduces toan effective-medium type approach. A DFT type dis-cussion is used in the initial theory of effective media tojustify the use of a form of the total energy which is thenparameterized using a variety of empirical models con-taining fit parameters. So systematic generalizations for T ∼ E F or higher T are not available although low- T fitted forms exist. Here we examine the more systematicfoundations which are, however, of little use in evaluatingmany of the current models which should be regarded asempirical fits to data bases.The total energy of the electron-ion system is writtenin EMA as: E tot = X i F { n i ( ~R i ) } + (1 / X i,j φ ( ~R i − ~R j ) (44)Here F { n i ( ~R i ) } is a function of the electron density, with n i ( ~R i ) being electron densities at atomic sites ~R i . Also, φ ( ~R i − ~R j ) is a pair-potential mainly used to model therepulsive interactions among the atoms. The atomic den-sities are associated with an embedding energy which istaken to be a function of the local average one-body elec-tron density ¯ n , following DFT. The embedding energy at T = 0 as a function of the electron density has beentabulated by various authors [56]. A superposition ap-proximation is introduced, and the total energy is finallyexpressed as a sum of terms for the isolated atom energiesplus the embedding energy in the homogeneous electrongas, and a gradient expansion in the electron density isused to allow for the density changes at atomic sites. Noattempt has been made to take advantage of ion-DFTusing an ion-correlation functional.However, this approach turns out to be inadequate inmany ways, yielding wrong elastic constants etc [54]. Theembedded atom model (EAM) was a numerical improve-ment on the EMA.The EAM attempts to merely use thestructure of the energy expression for fitting to parame-terized forms similar to it. Below we consider a modernform known as the so-called ‘glue model’ [57].The parameter-free first-principles NPA pair-potentialfor Aluminum is compared with one of the original force-matched Al potentials that uses some 40 fit parameters4and some 32 constraints. The ‘glue’ model uses the form: E tot = X i U { X j n ( r ij ) } + (1 / X i,j φ ( r ij ) (45)Here, instead of the electron gas term used in the EMA,or in the NPA as given by Eq. 40, a ‘glue’ function U ( n )is used, while φ ( r ) denotes the pair-potential. Hence, thisaluminum potential depends on specifying three func-tions, φ ( r ) , n ( r ) and U ( n ( r )) that contain fit parametersand constraints. Here again, simplification of the prob-lem that could be achieved using a DFT ion-correlationfunctional to include many-ion effects is not invoked. Theforces obtained from this potential (defined in terms ofparameters contained in polynomials or Pad´e forms orother fit functions) are matched (by adjusting the fit pa-rameters) to those from first principles calculations fora large variety of atomic configurations and physical sit-uations (solids, liquids, clusters, surfaces, defects etc.),including those at finite T . The Al potential of Ercolessiet al. has been fitted to liquid Al simulations at 1000Kand 2200K (i.e., T /E F = 0 . g ( r ) ofthe ionic structure appearing in the ion-correlation func-tional of the NPA, viz., Eq. 39. The needed g ( r ) is gen-erated in situ for fluids. In the case of specific crystalstructures at low T , one can use a known explicit form,as in Perrot’s Be calculations [25], or the phonon calcu-lations of Harbour et al [26].In Fig. 9 we compare the NPA Al-Al pair potentialwith that of the ‘glue potential’ of Ref. [57]. In panel(a) we display the NPA pair potential, labeled NPA2,for liquid aluminum at the normal density of 2.7 g/cm and T =2200K. The pair part of the glue potential is alsoshown (curve with triangles). Unlike in the case of theSW potential, the glue potential has Friedel like oscilla-tions, although not as long-ranged as in the NPA2 poten-tial. The force-matching methods do not usually have theaccuracy to recover the weaker Friedel oscillations whichdecay as 1 /r . The additional term of the glue potential, i.e, U ( n ) (not displayed in the figure) brings in contri-butions to the total energy that are also included in theNPA via the electron-fluid term and the XC-correlationenergies of ions as well as electrons, as indicated in Eq. 40.The ion-correlation energy is structure dependent, as in-dicated in Eq. 39.In Fig. 9 we also consider the case of Aluminum at adensity of 6.264 g/cm , and at a temperature of 1.75 eV.This corresponds to T /E F = 0 . 086 as ¯ Z = 3 for Alu-minum even at this temperature and compression. Thecompression drives up the Fermi energy and the effec-tive temperature T /E F remains low and the electronsare highly degenarate. The system has been studied ex-perimentally and theoretically by Fletcher et al. [58] us-ing a physically motivated ad hoc potential (YSRR6) [59]shown as a red dashed line. The YSRR potential is madeup of a Yukawa (Y) potential and a short-ranged repul-sive (SRR) potential. The YSRR6 potential at the den-sity of 6.24g/cm had been fitted to the DFT-MD g ( r ),and the corresponding S ( k ) is shown as a red-dashed linein panel (b). The label NPA6 refers to the first-principlespotential from the NPA at this density and temperature.It is a very different (black solid line) from the YSRRpotential, but gives excellent agreement with the MD-DFT S ( k ) as well as with the XRTS data of Fletcher. Itis found that the YSRR6 potential predicts a very lowcompressibility, very stiff phonon spectra, as well as un-realistic electrical conductivities, as discussed in detail byHarbour et al [60]. These emphasize the fact that fittingto a structural feature or even several is no guarantee ofhaving approximated the physically correct potential. VI. DISCUSSION The idea of using an ion correlation functional fortreating ion-ion many body effects was proposed as earlyas 1982 [27]. It was motivated by the success of using anelectron XC-functional to treat electron many-body ef-fects. The method which leads to a formulation in termof a pseudoatom, but it has not found significant adop-tion, partly because materials science has been mostlyconcerned with treating materials-engineering problemsinvolving defects, dislocations etc., or with complex ma-terials. Similarly, the idea of an ion correlation functionalis of not much use for molecules and other finite systemsused in quantum chemistry.However, there is little excuse for not using pseudo-atom approaches for uniform systems at finite- T as theconceptual and computational advantages are very sig-nificant. A common misconception, partly arising fromplacing excessive weight on the ‘chemical picture’ ofbonding leading to complex correlations, as well as thewish to subsume the electron subsystem in bonds, haslead to the belief that the NPA approaches do not con-tain, and cannot include, many-ion effects associatedwith the local structure of the medium. Here we havegiven examples showing how the NPA calculations in-5 r [Å] -0.100.1 V ii (r) [ e V ] NPA6YSRR6GLUENPA2 k [Å -1 ] S ii ( k ) a) b) NPA6: 6.26g/cm NPA2: 2.70g/cm FIG. 9. (Online color) Panel (a) displays the first-principlespair potential form NPA labeled NPA2 (blue dashed line) at T =2200K and the normal Al density of 2.7g/cm , as well asthe pair part of the ‘glue’ potential (triangles) of Ercolessi etal. NPA6 is the pair-potential for compressed Al at a den-sity of 6.26 g/cm and T =1.75 eV studied experimentally byFletcher et al [58]. A model potential (a Yukawa potentialadded to an ad hoc short-ranged repulsive potential) labeledYSRR6 used by Fletcher et al. is also shown (red dashedlines). In panel (b) the structure factors from the NPA6 po-tential and the YSRR6 potential are displayed. The YSRR6potential reproduces the experimental S ( k ) as it has been fit-ted to its g ( r ). The NPA6 is a first principles calculation withno fit parameters [60]. corporate such effects, and clarified the kinship of theNPA approach and its energy functonals to those used ineffective atom models. The NPA remains a strictly first-principles DFT approach, while the EMA approacheshave now become largely spring boards for data fitting.The ability of the NPA to study subtle effects like liquid-liquid phase transitions in complex materials like liquidSi and liquid C has been demonstrated and these havebeen reviewed in the context of how ion correlations areincluded in the NPA pair-potential scheme.Given that a very large data base of DFT results has been fitted the available EMA potentials, onemight envisage using them to construct first principlespair-potentials as well as a sophisticated first-principlesion-correlation functional, given that the electron XC-functional used in the primary DFT calculations is welldeveloped.Furthermore, although much effort has been directedto constructing kinetic energy functionals in the hope ofsimplifying DFT calculations, such calculations for N -atom systems will not yield useful ‘one-atom’ propertieslike ¯ Z , fractional compositions of ionic species, pseudopo-tentials and pair-potentials that are valuable interme-diates for additional calculations of physical properties.The NPA method can readily provide such quantities,and also easily include effects like van der Waals ener-gies, and quantum nuclear corrections in its potentials. VII. CONCLUSION The NPA method has been successfully applied to avariety of warm-dense-matter systems, as well as to anumber of solid-state systems. The detailed study of Sifor supercooled and hot liquid silicon has shown thatits accuracy is similar to DFT-MD simulations thatmay themselves differ from each other at the level ofthe electron XC-functionals used. The NPA has theadvantage that its simpler one-body electron densitiesallow the use of the local density approximation in theXC functional. The NPA yields one-atom propertieslike the ionization state, fraction compositions of ionicspecies, pseudopotentials etc., not easily available fromthe N-atom VASP-type simulations. On the other hand,the NPA being a static DFT theory, does not expose thecomplex bonding structure that may prevail transientlyin short-time scales. The thermodynamic properties ofthe system do not depend on such transient effects, andthe NPA deals only with their time average. DATA AVAILABILITY All the data used in this paper are available within thearticle in graphical form in the figures 1 to 9. If thereis any diffuclty in extracting them from the figures, thedata can be provided on request from the authors. [1] E. K. U. Gross, and R. M. Dreizler, Density FunctionalTheory , NATO ASI series, , 625 Plenum Press, NewYork (1993).[2] F. Perrot, and M.W.C. Dharma-wardana, Phys. Rev. E. , 5352 (1995).[3] M. W. C. Dharma-wardana, Phys. Rev. B , 155143(2019) DOI: 10.1103/PhysRevB.100.155143[4] M. W. C. Dharma-wardana, D. Neilson and F.M. Peeters Phys. Rev. B , 035303 (2019)https://link.aps.org/doi/10.1103/PhysRevB.99.035303DOI:10.1103/PhysRevB.99.035303[5] M. W. C. Dharma-wardana and F. Perrot Phys. Rev. Bvol , 035308 (2004). [6] J. N. Glosli Phys. Rev. Lett. , 4659 (1999).[7] M. W. C. Dharma-wardana, Contrib. Plasma Phys. , 135701 (2002)[9] M.S.Daw, S.M. Foiles and M.I. Baskes, , 251-310 (1993).[10] F.H. Stillinger and T. A Weber, Phys. Rev. B , 5262-5271 (1985).[11] J. Tersoff, Phys. Rev. B , 214103, (2005).[13] J. E. Saal, S. Kirklin, M. Aykol, B. Meredig, C. Wolver-ton, JOM , 1501 (2013). [14] S. Plimpton, Fast parallel algorithms for short-rangemolecular-dynamics, J. Comput. Phys. , 1-19 (1995).[15] S. Hamel, Lorin X. Benedict, Peter M. Celliers, M. A.Barrios, T. R. Boehly, G. W. Collins, Tilo D¨oppner, J.H. Eggert, D. R. Farley, D. G. Hicks, J. L. Kline, A.Lazicki, S. LePape, A. J. Mackinnon, J. D. Moody, H. F.Robey, Eric Schwegler, and Philip A. Sterne, Phys. Rev.B , 094113 (2012).[16] X.Gonze and C. Lee, Computer Phys. Commun. ,2582-2615 (2009).[17] G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996).[18] Valentin V. Karasiev, James Dufty, S.B. Trickey, Phys.Rev. Let. , 076401 (2018).[19] Olivier Hardouin Duparc, Philosophical Magazine Vol-ume , 3117-3131 (2009)[20] M.W.C. Dharma-wardana, Dennis D. Klug, and RichardC. Remsing Phys. Rev. Lett. , 075702 (2020) doi:10.1103/PhysRevLett.125.075702[21] M.W.C. Dharma-wardana in Laser Interactions withAtoms, Solids, and Plasmas , Carg’ese NATO workshop,1992, Edited by R.M. More (Plenum, New York, 1994),p311[22] L Harbour, and G. D. F¨orster, M. W. C. Dharma-wardana and Laurent J.Lewis, Physical review E ,043210 (2018).[23] John Ziman, Proc. Phys. Soc , 701 (1967)[24] L. Dagens, J. Phys. C: Solid State Physics, , 2333 (1972)[25] F. Perrot, Phys. Rev. E , 570 (1993).[26] L. Harbour, M. W. C. Dharmawardana, D. D. Klug, L.J. Lewis, Phys. Rev. E , 043201 (2017).[27] M. W. C. Dharma-wardana and F. Perrot, Phys. Rev. A , 2096 (1982)[28] C. E. Starrett, D. Saumon, J. Daligault, and S. HamelPhys. Rev. E , 033110 (2014).[29] M. S. Murillo, Jon Weisheit, Stephanie B. Hansen, andM. W. C. Dharma-wardana, Phys. Rev. E , 063113(2013).[30] B. Rosznayi, High Energy Density Physics 64, (2008)[31] P.A. Sterne S.B. Hansen, B.G. Wilson, W.A. Isaacs,HEDP, , 278 (2007)[32] F. Perrot and M. W. C. Dharma-wardana, Phys. Rev. B , 16536 (2000); Erratum: , 79901 (2003); arXive-1602.04734 (2016).[33] Richard C. Remsing, Michael L. Klein and Jianwei Sun,Physical Review B , 76 (1990).[35] M.W.C. Dharma-wardana and G.C. Aers, Phys. Rev. B. , 1701 (1983). [36] M. W. C. Dharma-wardana, Contr. Plasma Phys, , 85(2015).[37] F. Lado, S. M. Foiles and N. W. Ashcroft, Phys. Rev. A26 , 2374 (1983)[38] P. R. Smy, Advances in Physics, (5), 517-553, (1976)[39] S. H. Glenzer and Ronald Redmer, Rev. Mod. Phys. , 036407(2012).[41] M. W. C. Dharma-wardana, https://arxiv.org/abs/1705.10221 (2017)[42] James P. Edwards, Urs Gerber, Christian Schubert,Maria A. Trejo, Axel Weber Progress of Theoretical andExperimental Physics, (8), 083A01 (2017)[43] F. J. Rogers, H. C. Graboske Jr., and D. J. Harwood,Phys. Rev. A , 1577 (1970).[44] L. G. Stanton and M. S. Murillo, Phys. Rev. E , 043203(2016).[45] F. Perrot, Phys. Rev. A , 8334 (1991)[46] V. Recoules J. Clerouin, G. Zerah, P. M. Anglade, S.Mazevet, Phys. Rev. Lett. , 055503 (2006)[47] A. C. Maggs and N. W. Ashcroft, Phys. Rev. Lett. ,113 (1987).[48] K. Berland, V. R. Cooper, K. Lee, E. Schr‘¨odee, T. Thon-hauser, P. Hyldgaard, B. I. Lundqvist. Rep. Prog. Phys. (6), 066501 (2015)[49] Lucas J. Stanek, Raymond C. Clay III, M.W.C. Dharma-wardana, Mitchell A. Wood, Kristian R.C. Beckwith,Michael S. Murillo. https://arxiv.org/abs/2012.06451[50] J. Sun, B. Xiao, Y. Fang, R. Haunschild, P. Hao, A.Ruzsinszky, G. I. Csonka, G. E. Scuseria, and J. P.Perdew, Phys. Rev. Lett. , 106401 (2013).[51] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev.Lett. , 3865 (1996).[52] W. Yu, Z. Q. Wang, and D. Straoud, Phys. Rev. B ,13946 (1996)[53] F. Nadin, G. Jacucci and M.W.C. Dharma-wardana,Phys. Rev. A , 1025-1028 (1988). [NRCC 28307][54] Phys. Rev. B , 6443 (1994).[55] M. Finnis, J. Sinclair, Philosophical Magazine A ,583 (1994).[58] L. B. Fletcher et al , Nature Photonics , 274 (2015).[59] K. W¨unsch, J. Vorberger, and D. O. Gericke, Phys. Rev.E 79, 010201 (2009).[60] L. Harbour, M. W. C. Dharma-wardana, D. Klug and L.Lewis, Physical Review E94