Bypassing the malfunction junction in warm dense matter simulations
aa r X i v : . [ phy s i c s . c h e m - ph ] N ov Bypassing the malfunction junction in warm dense matter simulations
Attila Cangi and Aurora Pribram-Jones Max Planck Institute of Microstructure Physics, Weinberg 2, 06120 Halle (Saale), Germany ∗ Department of Chemistry, University of California, Irvine, CA 92697-2025, USA (Dated: September 11, 2018)Simulation of warm dense matter requires computational methods that capture both quantum andclassical behavior efficiently under high-temperature, high-density conditions. Currently, densityfunctional theory molecular dynamics is used to model electrons and ions, but this method’s com-putational cost skyrockets as temperatures and densities increase. We propose finite-temperaturepotential functional theory as an in-principle-exact alternative that suffers no such drawback. Wederive an orbital-free free energy approximation through a coupling-constant formalism. Our den-sity approximation and its associated free energy approximation demonstrate the method’s accuracyand efficiency.
PACS numbers: 71.15.Mb,31.15.E-
Warm dense matter (WDM) is a highly energeticphase of matter with characteristics of both solids andplasmas[1]. The high temperatures and pressures neces-sary for creation of WDM are present in the centers ofgiant planets and on the path to ignition of inertial con-finement fusion capsules[2, 3]. The high cost of experi-ments in this region of phase space has led to renewedinterest and great progress in its theoretical treatment[4–6]. Since both quantum and classical effects are crucialto accurate WDM simulations[7], density functional the-ory molecular dynamics has been used with increasingfrequency[8]. This method relies on Kohn-Sham densityfunctional theory (KS-DFT), which simplifies solving theinteracting many-body problem of interest by mappingit onto a non-interacting system[9, 10]. While the agree-ment between these calculations and experimental resultsis excellent[11, 12], the calculations are still incrediblyexpensive[13, 14]. The computational bottleneck in thesecalculations is the cost of solving the KS equations, a stepthat becomes increasingly expensive as temperatures riseand high-energy states become fractionally occupied. Infact, the computational cost exhibits nearly exponentialscaling with temperature due to the KS cycle including ahigh number of states at the temperatures of WDM[15].A solution to this problem is orbital-free DFT[16],which avoids this costly step by using non-interactingkinetic energy approximations that depend directly onthe electronic density. Because the kinetic energy is sucha large fraction of the total energy, however, these ap-proximations must be highly accurate to be of practi-cal use. Though much progress has recently been madefor WDM[17–19], approximations for thermal ensemblesare complicated by temperature effects. The KS ken-tropy, the free energy consisting of the non-interactingkinetic energy and entropy, must be approximated di-rectly, greatly complicating the production of useful and ∗ [email protected] D en s i t y xKS (0 K)KSPFATF FIG. 1. Shortcomings of the TF approximation in theWDM regime: Total density of five particles in the poten-tial v ( x ) = − ( πx/
10) within a box (of size 10 a.u.) atΛ = 1 / ( βµ ) = 0 .
93. Compare the exact density (solid blackcurve) with our PFA (dashed red curve) derived in Eq. (13),which is basically on top of the exact result. On the otherhand, the TF approximation (dotted green curve) capturesthe general qualitative features, but completely misses thequantum oscillations. We also show the corresponding exactdensity at zero temperature (light blue shaded area), with itspronounced oscillations that smooth as temperatures rise. efficient approximations.At zero temperature, potential functional theory(PFT) is a promising, alternative approach to the elec-tronic structure problem[20, 21]. It is also orbital-free,but skirts the troublesome issue of separately approxi-mating the KS kinetic energy. PFT’s coupling-constantformalism automatically generates a highly accurate ki-netic energy potential functional approximation (PFA)for any density PFA[20]. In this way, one needs only tofind a sufficiently accurate density approximation, as hasbeen demonstrated for simple model systems[22]. Ap-proximations to the non-interacting density have beenderived in various semiclassical[22–24] and stochasticapproaches[25]. It has also been shown that, at zerotemperature, PFT generates leading corrections to lo-cal approximations[22], which become exact in the well-known Lieb limit[26]. Finite-temperature Thomas-Fermitheory[27, 28] has been shown to become relatively exactfor non-zero temperatures under scaling similar to thatused by Lieb[29]. In this way, our method provides apathway to systematic improvements in approximationsto the kentropy, something generally missing from DFTapproaches.The particular scaling conditions under which TF be-comes exact for all temperatures is related to the break-down of purely quantum or purely classical behavior asboth temperatures and particle numbers increase[1]. Theimportance of both these effects in the WDM regime un-derlies its theoretical complexity[30], so it is useful torepresent the influences of temperature and density witha single parameter, an electron degeneracy parameter de-fined by Λ = 1 / ( βµ ), which depends on the system tem-perature 1 /β and temperature-dependent chemical po-tential µ . Then, the WDM regime can be defined aswhere Λ ≈
1. At these conditions, ground-state DFT ishugely expensive, while traditional plasma methods misscritical electronic structure features. In Fig. 1, we showthe density oscillations still present at WDM conditionsthat are neglected by the smooth, classical TF approxi-mation, but are captured by our method.In this work, we (i) derive PFT for thermal ensem-bles, (ii) give an explicit equation for the kentropy relyingsolely on the temperature-dependent density, (iii) deriveand implement a highly accurate density approximationin one dimension to illustrate our general result that con-tains the leading correction TF theory at finite tempera-ture, and (iv) perform (orbital-free) PFT calculations inthe WDM regime. Our method generates highly accuratedensity and kentropy approximations, skirts the need forseparate kentropy approximations, provides a roadmapfor systematically improved approximations, and con-verges more quickly as temperatures increase while main-taining accuracy at low temperatures. At the same time,it bridges low and high temperature methods, and so isuniquely suited to WDM.At non-zero temperature, the energy is replaced by thegrand canonical potential as the quantity of interest[31,32]. The grand canonical Hamiltonian is writtenˆΩ = ˆ H − β ˆ S − µ ˆ N , (1)where ˆ H , ˆ S , and ˆ N are the Hamiltonian, entropy, andparticle-number operators. In electronic structure the-ory, we typically deal with non-relativistic, interactingelectrons, most commonly within the Born-Oppenheimerapproximation. The electronic Hamiltonian (in atomic units here and thereafter) readsˆ H = ˆ T + ˆ V ee + ˆ V , (2)where ˆ T denotes the kinetic energy operator, ˆ V ee the in-terelectronic repulsion, and v ( r ) denotes the static exter-nal potential in which the electrons move. (We suppressthe spin of the electron for simplicity of notation.) Thegrand canonical potential can be written in terms of po-tential functionals as follows:Ω β v − µ = F β [ v ] + Z d r n β [ v ]( r )( v ( r ) − µ ) (3)with F β [ v ] = F β [ˆΓ v − µ ] = T [ˆΓ v − µ ]+ V ee [ˆΓ v − µ ] − β S [ˆΓ v − µ ]denoting the universal functional[33].In practice, approximating this direct expression wouldrequire two separate approximate potential functionals,one for the universal finite-temperature functional andone for the density:˘Ω β, dir v − µ = ˘ F β [ v ] + Z d r ˘ n β [ v ]( r )( v ( r ) − µ ) . (4)However, we can generate an approximation (denoted bya breve above the approximated quantity) to the univer-sal functional that corresponds to any density approx-imation. In analogy to the zero-temperature case[20],we introduce a coupling constant λ in the one-body po-tential, v λ ( r ) = (1 − λ ) v ( r ) + λv ( r ), where v is somereference potential. Via the Hellmann-Feynman theoremwe can then rewrite the grand potential,Ω β v − µ = Ω β + Z dλ Z d r n β [ v λ ]( r )∆ v ( r ) , (5)where ∆ v ( r ) = v ( r ) − v ( r ). Setting v = 0 and defining¯ n β [ v ]( r ) = R dλ n β [ v λ ]( r ), we can now write the finite-temperature universal functional in terms of the densitywritten as a potential functional: F β, cc n β [ v ] = Z d r { ¯ n β [ v ]( r ) − n β [ v ]( r ) } v ( r ) . (6)This defines an approximate functional, ˘ F β, cc ˘ n β [ v ], corre-sponding to the chosen density approximation ˘ n β and isthe generalization of PFT to thermal ensembles.Since practical use of this formula as written wouldrequire sufficiently accurate approximations to the inter-acting electron density (which are likely unavailable), weinstead apply it to the non-interacting electrons of theKS system. In DFT, the KS system is a clever way ofapproximating the exact F β by mapping the interact-ing system to a non-interacting system with the sameelectronic density and temperature. This determinesthe one-body KS potential and corresponding chemicalpotential. Through this mapping, the non-interacting,finite-temperature universal density functional can bedefined[33] F β s [ n ] := min ˆΓ → n K β [ˆΓ] = K β [ˆΓ β s [ n ]] = K β s [ n ] , (7)which generates the KS equations and, through their so-lutions, the KS orbitals. Slater determinants of the KSsingle-particle orbitals are the KS wavefunctions. Theorbitals are implicit functionals of the density via the KSequations, and the average density can be constructedby Fermi-weighted summing of the orbitals. Solution ofthese equations at every time-step is the most costly stepof DFT molecular dynamics.The KS potential is defined[20, 21] v S ( r ) = v ( r ) + ˜ v H [ n β S [ v S ]]( r ) + ˜ v XC [ n β S [ v S ]]( r ) , (8)where, in contrast to KS-DFT, the density is posed as a potential functional, and tildes distinguish density func-tionals from potential functionals. All many-body inter-actions among the electrons are captured in the usual KS-DFT sense, via the (traditionally defined) Hartree andXC potentials[34]. The difference from a usual KS-DFTcalculation is that Eq. (8) in conjunction with an ap-proximation to the non-interacting density bypasses thehugely expensive iterative solution of the KS equationsfor WDM. Choosing a potential functional approximationto the non-interacting density automatically generates anapproximated KS potential, as illustrated in the Supple-mental Materials. Once the self-consistent KS potentialis determined, the KS kentropy is computed from K β, ccS ,n β S [ v S ] = Z d r { ¯ n β S ( r ) − n β S [ v S ]( r ) } v S ( r ) , (9)which is the analog of Eq. (6) for KS electrons.Again, Eq. (9) defines a coupling-constant approxima-tion, ˘ K β, ccS , ˘ n β S [ v S ], when evaluated on any chosen approx-imation to the non-interacting density ˘ n β S . Finally, thegrand potential expressed in terms of KS quantities[33],Ω β v − µ = K β, ccS ,n β S [ v S ] + U [ n β S [ v S ]] + F β xc [ n β S [ v S ]]+ Z d r n β [ v S ]( r ) ( v ( r ) − µ ) , (10)can be evaluated via Eq. (9). Through this result, weleverage the body of time-proven XC approximations andeliminate the need to construct separate approximationsto the KS kentropy for use in orbital-free (and therebycomputationally inexpensive) schemes. Only an approx-imation to the non-interacting density is required. Thismeans that a general, systematic, non-empirical route toimproved kentropy approximations is now available.In principle, a possible starting point for deriving anapproximation to the non-interacting density at finitetemperature is the semiclassical propagator G β sc ( r , r ′ ; α ) = G ( r , r ′ ; α ) f β ( α ) , (11) which can be written as a convolution of the zero-temperature propagator G ( r , r ′ ; α ) with the factor f β ( α ) = πα/ [ β sin( πα/β )] carrying all temperature de-pendence, and α is a complex variable[35]. From thepropagator, we extract the density via an inverse Laplacetransformation[36]˘ n β S [ v S ]( r ) = lim r ′ → r πi η + ∞ Z η −∞ dα e µα α G β sc [ v S ]( r , r ′ ; α ) . (12)To illustrate the significance of our main result inEq. (9), we consider a simple, yet useful, numericaldemonstration: Non-interacting, spinless fermions in anarbitrary potential v ( x ) confined to a box of size L obey-ing vanishing Dirichlet boundary conditions. (In a prac-tical realization, this would be the self-consistent KS po-tential of the given many-body problem.) Recently, ahighly accurate PFA to the density was derived for thismodel using the path integral formalism and semiclassi-cal techniques[37]. Here we extend this result to finitetemperature via Eq. (11) and obtain:˘ n β S ( x ) = lim x ′ → x X λ =1 ∞ X j =0 ˘ γ β S ( x, x ′ ; λ, j ) , (13)a PFA to the density at a given temperature and chemicalpotential, where˘ γ β S ( x, x ′ ; λ, j ) = sin Θ λµ ( x, x ′ ; j )csch[ πβ T λµ ( x, x ′ ; j )]( − λ +1 β p k µ ( x ) k µ ( x ′ ) . (14)Here we define generalized classical phases Θ µ ( x, x ′ ; j ) = θ − µ ( x, x ′ ) + 2 jθ µ ( L ), Θ µ ( x, x ′ ; j ) = θ + µ ( x, x ′ ) + 2 jθ µ ( L ),Θ µ ( x, x ′ ; j ) = θ − µ ( x, x ′ ) − j + 1) θ µ ( L ), Θ µ ( x, x ′ ; j ) = θ + µ ( x, x ′ ) − j + 1) θ µ ( L ) and generalized classical travel-ing times T λµ ( x, x ′ ; j ) = d Θ λµ ( x, x ′ ; j ) /dµ . Furthermore, θ ± ( x, x ′ ) = θ ( x ) ± θ ( x ′ ), where θ µ ( x ) = R x dy k µ ( y ) and k µ ( x ) = p µ − v ( x )) at a given chemical potential µ ,which is determined by normalization of the density.The physical interpretation of our result in Eq. (13)is instructive: For a given chemical potential there areinfinitely many classical paths that contribute to the to-tal density. The paths are classified into four primitives(identified by λ ) onto which an integral number of pe-riods (labelled by j ) is added. The first primitive isspecial, in the sense that it yields the TF density. Allother primitives and additional periods carry phase infor-mation about reflections from the boundaries, producingquantum density oscillations that greatly improve uponthe TF result[37]. For a more details on the derivationand the physical interpretation, we refer to Ref. [37].Our result in Eq. (13) can be evaluated numericallyfor a given temperature by truncating the infinite sumat a conveniently chosen upper limit such that the re-sult is convergent. Importantly for WDM applications,the higher the temperature, the lower the upper limit re-quired for convergence of the sum. In fact, in the WDMregime only the leading term ( j = 1) in the sum needs tobe kept. While especially powerful at finite temperature,this might be a universal feature due to the semiclassicalnature of our approximation. Similar results have alsobeen recently found at zero temperature[37, 38].However, the stationary phase approximationused to derive Eq. (13) yields the TF densityat zero temperature as the leading term, i.e.,lim x ′ → x ˘ γ β S ( x, x ′ ; 1 ,
0) = k µ ( x ) /π = ˘ n TF ( x ), in-stead of the finite-temperature TF density ˘ n β TF ( x ) = p / (2 π β ) F ( z ), where F ( z ) = R ∞ da {√ a [1 + exp( a − z )] } − and z = β k µ ( x ) /
2. We fix this problem withan ad-hoc correction and ensure the correct boundaryconditions. To do so, we replace the density from thefirst primitive lim x ′ → x ˘ γ β S ( x, x ′ ; 1 ,
0) with a Gaussianinterpolation of ˘ n TF ( x ) and ˘ n β TF ( x ). In this way, wecope with the density approaching the high-temperaturelimit (under which TF theory becomes exact) differentlyin two distinct regions, the interior of the box and theedge regions near the walls. These may be considered astwo distinct boundary layers with different asymptoticexpansions in the high-temperature limit. Note that thesize of the boundary layers in the edge regions shrinksas the limit is approached. The Gaussian interpola-tion applied here is a crude version of the approachused in boundary-layer theory to match two differentapproximations with different asymptotic behavior[39].In Fig. 1, we plot a typical density of five particlesin the WDM regime (Λ ≈
1) in the potential v ( x ) = − ( πx/
10) within a ten-unit box, along with approx-imate densities. The black curve is the exact result, thered dashed curve is our approximation, and the greendotted curve is the TF density. In addition, the light-blueshaded area denotes the corresponding density at zerotemperature. This figure demonstrates that quantum os-cillations in the density persist in the WDM regime andthat TF theory completely fails to capture them. On theother hand, our PFA – derived to include quantum ef-fects – is able to describe them properly and is thereforehighly accurate.Next, we demonstrate the accuracy of our approach forkentropies. For our example, Eq. (9) simplifies to˘ K β, ccS , ˘ n β S [ v ] = K β S , + Z dx (cid:8) ˘¯ n β S ( x ) − ˘ n β S [ v ]( x ) (cid:9) v ( x ) . (15)In this case the reference potential is not zero, but an infi-nite square well. Hence, a kentropic contribution K β S , = T β S , − S β S , /β of the reference system appears, which wecan compute exactly quite simply. The kinetic energyof the infinite square well is T β S , = P Nj f β j ǫ j, , and theentropy is S β S , = − P j f β j ln( f β j ) + (1 − f β j ) ln(1 − f β j ),with f β j = 1 / (1 . β ( ǫ j, − µ )]) denoting Fermifunctions and ǫ j, and µ the j th eigenvalue and chem-ical potential for the non-interacting reference system. TABLE I. Residual kentropy of five particles in the same po-tential as in Fig. 1. We list the error of the conventional TFapproach and of our PFA (given in Eq. (15)) far below andabove where WDM is typically encountered.Λ K β S , ∆ K β S error × TF PFA0 .
16 3 .
94 0 .
462 6 . − . .
31 3 .
87 0 .
461 7 . − . .
47 3 .
76 0 .
459 7 . − . .
62 3 .
64 0 .
456 8 . − . .
78 3 .
50 0 .
452 8 . − . .
93 3 .
34 0 .
448 8 . − . .
09 3 .
16 0 .
444 8 . − . .
40 2 .
77 0 .
435 8 . − . .
71 2 .
36 0 .
425 7 . − . .
02 1 .
92 0 .
414 7 . − . .
48 1 .
25 0 .
396 6 . − . .
94 0 .
58 0 .
378 5 . − . . − .
10 0 .
360 5 . − . . − .
99 0 .
338 4 . − . We avoid temperature-dependent KS eigenvalues[30] bychoosing a purely non-interacting reference system, not aKS system associated with a specific interacting system.Evaluating Eq. (15) for the same potential as in Fig. 1yields the results in Tab. I. We measure the error of TFtheory and our PFA with respect to the residual kentropy∆ K β S = K β, ccS ,n β S − K β S , , because this is the only piece of thetotal kentropy being approximated. We also list the elec-tron degeneracy parameter values and the kentropic con-tributions of the reference system. Reaching from coldtemperatures up to the WDM regime (around Λ ≈ -0.05 0 0.05 0.1 0.15 0.2 0 2 4 6 8 10 R e s i dua l k en t r op i c den s i t y xKSPFATF FIG. 2. Residual kentropic density of five particles in the samepotential as in Fig. 1 in the WDM regime. Our PFA (solidred curve) derived in Eq. (13) is on top of the exact result(solid black curve). The TF result (dotted green curve), onthe other hand, follows the general trend as expected, butmisses quantitative details. for WDM, where solution of the KS equations for nu-merous occupied states becomes especially daunting. Weretain the advantages of the KS system while avoidingthe costly, repetitive solution of eigenvalue problems byisolating a small piece of the kentropy to approximatethrough the coupling-constant formalism. The referencesystem is always chosen such that its kentropy is knownexactly. Combined with our density approximation, thisimproves approximate kentropies by up to an order ofmagnitude in the WDM regime and produces highly ac-curate kentropic densities. The density approximationderived in this paper is computationally efficient becauseonly the leading terms are needed for convergence atWDM temperatures.The path integral method used to derive thisapproximation[37] invites use of successful zero-temperature approximations to the propagator, and it isa promising approach for extension to three dimensionalsystems. Furthermore, combining the presented finite-temperature PFT with semiclassical methods offersprospects for a systematic route to exchange energyapproximations, instead of only relying on existing,zero-temperature density functional approximations.Work in this direction is currently in development. Withthese advantages, finite-temperature PFT is poised tobridge the “malfunction junction” of WDM by providingcomputationally efficient, semiclassical methods at hightemperatures and densities.
ACKNOWLEDGMENTS
We acknowledge Hardy Gross and Kieron Burke forproviding a fruitful atmosphere facilitating independentresearch. We are grateful to Rudy Magyar for usefuldiscussion. A.C. has been partially supported by NSFgrant CHE-1112442. A.P.J. is supported by DOE grantDE-FG02-97ER25308. [1] F. Graziani, M. P. Desjarlais, R. Redmer, and S. B.Trickey, eds.,
Frontiers and Challenges in Warm DenseMatter , Lecture Notes in Computational Science andEngineering, Vol. 96 (Springer International Publishing,2014).[2] S. Atzeni and J. Meyer-ter Vehn,
The Physics of InertialFusion: Beam-Plasma Interaction, Hydrodynamics, HotDense Matter (Clarendon Press, 2004).[3] N. R. C. C. on High Energy Density Plasma PhysicsPlasma Science Committee,
Frontiers in High EnergyDensity Physics: The X-Games of Contemporary Science (The National Academies Press, 2003).[4] T. R. Mattsson and M. P. Desjarlais, Phys. Rev. Lett. , 017801 (2006).[5] F. R. Graziani, V. S. Batista, L. X. Benedict, J. I. Castor,H. Chen, S. N. Chen, C. A. Fichtl, J. N. Glosli, P. E.Grabowski, A. T. Graf, S. P. Hau-Riege, A. U. Hazi, S. A.Khairallah, L. Krauss, A. B. Langdon, R. A. London,A. Markmann, M. S. Murillo, D. F. Richards, H. A. Scott,R. Shepherd, L. G. Stanton, F. H. Streitz, M. P. Surh,J. C. Weisheit, and H. D. Whitley, High Energy DensityPhysics , 105 (2012).[6] K. Y. Sanbonmatsu, L. E. Thode, H. X. Vu, and M. S.Murillo, J. Phys. IV France , Pr5 (2000).[7] M. D. Knudson and M. P. Desjarlais, Phys. Rev. Lett. , 225501 (2009).[8] B. Holst, R. Redmer, and M. P. Desjarlais, Phys. Rev.B , 184201 (2008).[9] P. Hohenberg and W. Kohn,Phys. Rev. , B864 (1964).[10] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[11] A. Kietzmann, R. Redmer, M. P. Desjarlais, and T. R.Mattsson, Phys. Rev. Lett. , 070401 (2008).[12] S. Root, R. J. Magyar, J. H. Carpenter, D. L. Hanson,and T. R. Mattsson, Phys. Rev. Lett. , 085501 (2010).[13] D. Marx and J. Hutter, in Modern Methods and Algo-rithms of Quantum Chemistry , NIC Series, Vol. 1, editedby J. Grotendorst (J ulich, 2000) pp. 301–449.[14] D. Marx and J. Hutter,
Ab Initio Molecular Dynamics:Basic Theory and Advanced Methods (Cambridge Uni-versity Press, 2009).[15] V. V. Karasiev, T. Sjostrom, D. Chakraborty, J. W.Dufty, K. Runge, F. E. Harris, and S. Trickey, in
Frontiers and Challenges in Warm Dense Matter , Lec-ture Notes in Computational Science and Engineering,Vol. 96, edited by F. Graziani, M. P. Desjarlais, R. Red-mer, and S. B. Trickey (Springer International Publish-ing, 2014) pp. 61–85.[16] Y. A. Wang and E. A. Carter, in