Phase Diagram and Electronic Structure of Praseodymium and Plutonium

We develop a new implementation of the Gutzwiller approximation in combination with the local density approximation, which enables us to study complex 4f and 5f systems beyond the reach of previous approaches. We calculate from first principles the zero-temperature phase diagram and electronic structure of Pr and Pu, finding good agreement with the experiments. Our study of Pr indicates that its pressure-induced volume-collapse transition would not occur without change of lattice structure—contrarily to Ce. Our study of Pu shows that the most important effect originating the differentiation between the equilibrium densities of its allotropes is the competition between the Peierls effect and the Madelung interaction and not the dependence of the electron correlations on the lattice structure.


I. INTRODUCTION
There has been a renewed interest in first-principles approaches to the electronic structure of strongly correlated materials.Density functional theory (DFT)-and, in particular, the local density approximation (LDA)-proved to be a good starting point for deriving model Hamiltonians [1,2] that can be studied with more elaborate methods that are able to treat correlations.While early approaches to solve the realistic many-body problem in solids focused on perturbative treatments of the interactions [3], over the last two decades several nonperturbative methodologies have emerged.Dynamical mean field theory (DMFT) was combined with realistic electronic structure methods, for example, in the LDA þ DMFT approach [4,5].This methodology can be thought of as a spectral density functional [6] and is nowadays widely used to study 3d, 4d, 5d, 4f, and 5f systems.For reviews see, e.g., Refs.[7][8][9][10][11].LDA þ DMFT has been implemented in different basis sets, such as the linearized augmented plane wave (LAPW) [12,13], plane-wave pseudopotentials [14], the projector-augmented wave method [15], and linearized muffin-tin orbitals [16].Another important approachwhich is not as accurate as DMFT but has the advantage to be less computationally demanding-is the Gutzwiller approximation (GA) [17][18][19], which was first implemented to study real solids in Ref. [20].The GA approximation was, thereafter, extensively developed [21][22][23][24][25][26][27], and it has been formulated and implemented in combination with realistic electronic structure calculations such as the LDA þ GA approach [23,28], which has been applied successfully to many systems [29][30][31][32][33][34][35][36].A third important many-body technique is the slave boson approach (SB) [37,38], which is, in principle, an exact reformulation of the quantum many-body problem for model Hamiltonians, and it reproduces the results of the GA at the saddle-point level [24,39].This technique has recently been extended to treat full rotationally invariant interactions [38,40], and it has also been combined with LDA for the study of real materials, either in the form of impurity solvers for the resulting LDA þ DMFT impurity models (LDA þ DMFTþ SB) [6,41] or directly on the lattice [42]-which are equivalent approaches, as we will show.
On the methodological side, here we show that the three above-mentioned methods are closely connected and largely complementary of each other.We use the connection between the GA and the SB methods to introduce a functional formulation, which can be used not only at zero temperature but also at finite temperatures.This functional is a first step toward deriving formulas for the forces [43] and the phonon spectra in the LDA þ GA and LDA þ SB methods.Our functional formulation of the LDA þ GA method [23] has a mathematical structure similar to LDA þ DMFT [7].This parallelism suggests possible synergistic combinations between the two methods-such as using LDA þ GA for structural relaxation while using the exact impurity solver in the LDA þ DMFT iteration to determine the spectral properties.Furthermore, it enables us to pattern the LAPW interface [44] between LDA and the GA after the LDA þ DMFT work of Ref. [12].These connections result in a new algorithm for solving the LDA þ GA equations, which is faster and more precise than earlier methods and sheds light on the physical interpretation of the SB amplitudes-which are central quantities both in the GA and in the SB approach.In fact, we display a connection between the SB amplitudes and the coefficients of the Schmidt decomposition [45]-which was also recently used to derive the density matrix embedding theory [46][47][48][49].Our algorithm consists in recursively calculating the ground state of a series of Anderson impurity Hamiltonians (one for each inequivalent impurity within the lattice unit cell), whose baths have the same dimension as the corresponding impurities.This enables us to derive accurate equations of state for materials currently far beyond the reach of LDA þ DMFT.
The technical advances obtained in this work result in a new understanding of the volume-collapse transition in f systems.In particular, we use our all-electron (LAPW) implementation of the LDA þ GA method to study two prototypical systems with partially delocalized f electrons: elemental Pr and Pu.
Pr is a rare earth like Ce, and it is the next element in the periodic table.An interesting property of Pr is that, similarly to many other rare-earth compounds, it undergoes a volume-collapse structure transition under pressure, which is accompanied by an abrupt delocalization of the f electrons [50].Here, we compute its pressure-volume phase diagram, finding very good agreement with the experiments.In particular, we show that the method is able to capture the pressure-induced volume-collapse structure transition toward the low-symmetry α-U phase, and that the GA correction to the total energy is crucial to correctly determine the stable lattice structure of Pr.Finally, we investigate the relation between the f delocalization and the volume-collapse structure transition-which is one of the most important puzzles in condensed matter physics.Our main conclusion is that, contrarily to Ce [35], in Pr there would not be any volume-collapse transition without taking into account the change of lattice structure (at least at low temperatures).
The stable allotrope of Pu at ambient conditions is α-Pu, and five different crystalline phases (named β; γ; δ; δ 0 ; ϵ) can be stabilized at higher temperatures.One of the most intriguing properties of Pu is that these temperature-induced structure transitions are accompanied by significant changes of density.Here, we perform LDA þ GA calculations of all of the six phases of Pu and study how the total energy and the f-electron correlations depend on the volume and the lattice structure.These results provide a complete bird's eye view of this material and, in particular, indicate that the most important effect originating the above-mentioned large differences between the equilibrium volumes of the phases of Pu is the competition between the Peierls effect and the Madelung interaction and not the dependence of the electron correlations on the lattice structure, which we find to be a negligible effect.We point out that the explanation of this phenomenon is of great interest both physically and from a metallurgic standpoint.
The outline of the paper is as follows.In Sec.II, the formulation and the implementation of the GA/SB developed in Refs.[22,24,25,27] are substantially improved.In Sec.III, the connection between the GA/SB and DMFT is discussed.In Sec.IV, our functional formulation of the LDA þ GA method is derived.Finally, in Secs.V and VI, our calculations of Pr and Pu are illustrated.

II. THE GUTZWILLER APPROXIMATION FOR THE HUBBARD MODEL
Let us consider the Hubbard model (HM) where k is the momentum conjugate to the unit-cell label R, the atoms within the unit cell are labeled by i; j, and the spin orbitals are labeled by α; β.We assume that the first term is nonlocal, i.e., that and that Ĥloc includes both the one-body and the two-body part of the Hamiltonian.Note that no specific assumption needs to be made on the structure of the local interaction.
In particular, in this work we have used the rotationally invariant Slater-Condon parametrization of the on-site interaction [51].
We define the temperature T and the corresponding fermionic Matsubara frequencies as ω ¼ ð2m þ 1ÞπT .As shown in Appendix A, the SB mean-field theory [37,38] can be formulated in terms of the following Lagrange function, which gives the free energy when it is evaluated at its saddle point and reduces to the GA at T ¼ 0: where the dimension M i of the matrices R i ; D i ; λ i ; λ c i ; η i is the number of spin orbitals within the ðR; iÞ space; R, λ and η are block matrices whose blocks are R i , λ i and η i , respectively; μ is the chemical potential; N is the total number of electrons (normalized to the number of k points N ); the matrices F ia and H loc i represent the local operators fRia and Ĥloc Ri in a given (arbitrary) basis set of local multiplets jΓ; Rii, and ϕ i are the SB amplitudes, which are i.e., that couple only local multiplets with the same number of electrons (which amounts to restricting ourselves to nonsuperconducting phases).Note that the above finite-temperature extension has been recently used in Ref. [36] to study the thermodynamical properties of Ce in relation to its α-γ isostructural volume-collapse transition.
For later convenience, we observe that the first term of Eq. ( 3) can be equivalently rewritten as where and in the last step of Eq. ( 7), we used that the renormalization factors in the quasiparticle Green's function give only a frequency-independent term, which does not contribute because of the Matsubara summation.Note that since R, λ, and η are block matrices, ΣðzÞ is also block diagonal, with blocks Once the GA solution is determined by imposing the saddle-point conditions of L N with respect to all of its arguments, the expectation value of any observable can be readily computed.In particular (see Appendix A), it can be proven that the expectation value of any local operator within the site ðR; iÞ is given by where A i is the representation of ÂRi in the same basis used in Eqs. ( 4) and (5), A. Physical interpretation of the parameters ϕ i based on the Schmidt decomposition In this section, we show that the space of matrices ϕ i [see the second line of Eq. ( 3)] can be conveniently mapped into the Hilbert space of states jΦ i i of an impurity system composed of the i impurity and an uncorrelated bath with the same dimension, which is determined self-consistently in order to describe the entanglement between Ĥloc i and the rest of the system [see Eq. ( 1)].
Let us consider the impurity local many-body space V i generated by the same Fock fermionic basis that has been used in Eqs. ( 4) and (5), where ν Γ α are the occupation numbers of jΓ; ii, so that the total number of electrons corresponding to this state is We define a copy of W i of V i generated by another set of Fock states jn; ii represented as Note that, since the f ladder operators act only on the W i degrees of freedom, they anticommute with all of the c ladder operators.We define ν n α as the occupation numbers of jn; ii, so that the total number of electrons corresponding to this state is Let us now consider a generic pure state within the tensor-product space E i ≡ V i ⊗ W i , and represent it as where U PH is the particle-hole (PH) transformation satisfying the following identities: i.e., acting only on the f degrees of freedom.
For later convenience, we assume that the condition [Eq.(6)] is respected by the matrix ϕ i appearing in Eq. ( 16) so that it couples only states with N Γ ¼ N n .Note that this condition amounts to assuming that where is the total number operator in the embedding system E i , and M i is the number of spin orbitals in the R; i space.In fact [see Eq. ( 16)], jΦ i i is a linear combination of product states with It can be readily verified that, for any operator Â acting within the space generated by the states jΓ; ii, where ϕ i is the matrix of coefficients that appears in Eq. ( 16) and Furthermore, it can be shown that where the matrix elements of F iα are defined by and can be equivalently calculated as Note, in fact, that Eqs. ( 28) and ( 29) represent the matrix elements of ladder operators in their own Fock basis [see the definitions ( 12) and ( 14)].The explicit derivation of Eqs. ( 23)-( 27) is given in Appendix B. Thanks to Eqs. ( 23)-( 27), the GA Lagrange function [see Eq. ( 3)] can be rewritten as follows: where and jΦ i i belongs (by construction) to the subspace defined by Eq. ( 21); i.e., it is an eigenstate of Ntot i with eigenvalue M i .
Since the GA solution is stationary with respect to jΦ i i and E c i , Ĥemb i can be interpreted as an impurity Hamiltonian whose bath has the same dimension of the impurity and is determined by the GA procedure in order to describe the entanglement between the impurity and the rest of the system.
We point out that, in principle, a Hamiltonian whose bath has the same dimension of the impurity is sufficient to represent exactly the ground-state local properties of the Hubbard model, as it can be readily shown by making use of the Schmidt decomposition [45,46].What we have shown in this section is that the GA consists in assuming that the form of this Hamiltonian is limited to an Anderson impurity model [see Eq. ( 31)], where the coefficients D i ; λ c i are determined by the stationarity of the Lagrange function [Eq.(30)].This insight makes it clear that taking into account all of the components of ϕ i is crucial to accurately describe the local physics of the system-including the offdiagonal matrix elements in a basis that diagonalizes Ĥloc i .

B. Numerical solution of the GA Lagrange equations
In Eq. (7) and the text below, we have introduced the block matrices R, λ, and η appearing in the functional [Eq.(30)], whose respective blocks (one for each i within the unit cell) are R i , λ i , and η i .For later convenience, we define Π i as the projectors onto the above-mentioned i local subspaces.The symbol f will indicate the Fermi function.
It can be readily shown that the saddle-point condition of L N with respect to all of its arguments provides the following system of Lagrange equations: Note that the projectors Π i appear in Eqs. ( 32)- (35) because derivatives are taken with respect to the matrix elements of the block matrices η, λ i , and R i .
A possible way to compute the Gutzwiller solution is the following [27]: (i) Given ðR; λÞ, use Eqs.(32) and (33) to compute the Lagrange multipliers μ and η and the corresponding jΨ 0 i, which determines n 0 i through Eq. ( 34), D i through Eq. (35), and λ c i through Eq. ( 36).(ii) Thereafter, build the embedding Hamiltonians Ĥemb i and compute jΦ i i [see Eq. ( 37)], which determine the left sides of Eqs.(38) and (39).Equations (38) and (39) are verified if and only if ðR; λÞ is the correct set of variational parameters.
In conclusion, we have formulated the solution of the Gutzwiller equations as a root problem for ðR; λÞ, which can be formally written as and solved numerically.Note that the vector functions F i [see Eqs.(38) and (39)] can be evaluated independently through the numerical steps outlined above.

C. Summary of the main results of Sec. II
In this section, we have expressed the GA Lagrange function derived in Ref. [27] [see Eq. ( 3)] in the convenient form of Eq. (30), from which follows an exceptionally efficient numerical scheme [see Eqs. ( 32)-( 39) and text below].
We have shown that our algorithm consists in iteratively solving a series of Anderson impurity Hamiltonians whose bath has the same dimension as the impurity [see Eq. (31)].This finding provides a useful interpretation of the Gutzwiller variational parameters based on the Schmidt decomposition and opens up the possibility to exploit techniques such as those developed in quantum chemistry to solve Ĥemb i , in order to further speed up our algorithm.
We point out that in our numerical scheme, the treatment of the correlation effects scales linearly with the number of correlated atoms per unit cell (as in LDA þ DMFT).Since linear-scaling DFT methods are also available [52], the linear-scaling property of our solver opens up the possibility of studying correlated systems with extremely large supercells [53][54][55]-even containing several hundreds of correlated atoms.
In Appendix C, alternative expressions for the GA Lagrange equations ( 32)- (35) are derived using the Green's function formalism [see Eqs.(C8) and (C9)], which can be preferable if the unit cell of the system contains many atoms and/or if only a few orbitals are correlated [56].In Appendix D, the numerical strategy discussed in this section is generalized to Anderson impurity models.

III. THE GUTZWILLER-BAYM-KADANOFF FUNCTIONAL
In this section, we formulate the GA of the Hubbard model [see Eq. ( 1)] as the saddle-point of a functional of the coherent part of the local Green's function, and we show that the mathematical structure of this functional resembles the Baym-Kadanoff (BK) theory [57] on top of the DMFT approximation.
Let us rewrite Eq. ( 30) as follows: where we have implicitly assumed that the regularization factor e iω0 þ is present in the Matsubara summation, and Θ Ĥloc i depends on Ĥloc i through Eq. (31).Note that, thanks to the formal manipulations of Eq. ( 7), the quasiparticle parameters R; λ; η affect the first term of L 0 only through the Gutzwiller self-energy Σ, which was defined in Eq. (8).
It is useful to promote the self-energy to an independent variable by introducing the following additional Lagrange-Legendre term in Eq. ( 41): where G i ðiωÞ are, at the present stage, the Lagrange multipliers used to enforce the GA definition of Σ i ðiωÞ [see Eq. ( 9)].However, by deriving the so-obtained Lagrange function with respect to Σ i ðiωÞ, one obtains the Dyson equation for the i local Green's function, where Π i is the projector onto the i local subspace.
The above formal manipulations enable us to express the GA in terms of the following Lagrange function: where We observe that the functional Φ Ĥloc i ;N can be viewed as a Lagrange function on its own, which depends parametrically on G i and μ, and explicitly on The stationary solution of Φ Ĥloc i ;N for the variables X i can be formally expressed as a function of G i and μ itself.This mathematical construction enables us to define the following functional of G i and μ only: which can be substituted back in Eq. (45).
In conclusion, we have demonstrated that the GA solution is the saddle-point with respect to G and Σ of the functional, which resembles the BK theory on top of the DMFT approximation.
We point out that, remarkably, the functional Φ GA depends on the nonlocal dispersion ϵ k only through the coherent part of the local Green's functions G i .In other words, it is determined only by Ĥloc i and N, and is formally the impurity Luttinger functional corresponding to the i onsite interaction, and T P ω Tr½G i ϵ loc i is the additional term arising from having included the on-site quadratic part ϵ loc of the Hamiltonian and the chemical potential μ within the definition of the self-energy.Note that, since the on-site quadratic operators have to be treated together with the interaction within the GA, the Gutzwiller approximation for Φ L Ĥint i alone cannot be defined in general.

A. Summary of the main results of Sec. III
In this section, we have derived a functional formulation of the GA which has essentially the same formal structure of the Baym-Kadanoff theory on top of the DMFT approximation [see Eq. ( 48)], making an exception for the following technical differences.(1) The Gutzwiller-Baym-Kadanoff functional depends only on the coherent part of the Green's function [see Eq. ( 44) and Eqs.(C5)-(C7)].(2) Within the GA, it is necessary to treat the quadratic part of the local Hamiltonian together with the interaction [see Eq. ( 2) and text below Eq. ( 48)].This result clarifies the connection between GA/SB and DMFT and proves the equivalence between DMFT þ SB and SB.

IV. FUNCTIONAL FORMULATION
OF LDA þ GA Approximations to DFT [58] represent the state of the art of materials simulations.DFT calculations based on LDA [59] enable us to theoretically attack a wide class of materials, but they are generally not satisfactory for the socalled "strongly correlated" systems, such as, e.g., the high-T c superconductors, transition metal oxides, and rare-earth compounds.In order to study this important class of materials, several "hybrid" techniques, such as LDA þ U [51], LDA þ DMFT [7], and LDA þ GA [23,28], have been developed.
In this section, we discuss our implementation of the LDA þ GA method.

A. Correlated orbitals
The application of any LDA þ X method requires the identification of a proper subset of "correlated" orbitals (e.g., d or f), which we indicate with the symbol P.
The correlated orbitals are determined on a physical basis.Usually, they are constructed using Wannier function methods, e.g., maximally localized Wannier functions [60], projected Wannier functions [16,61], or quasi-atomic minimal basis-set orbitals [62,63].In particular, in this work, we refer to the construction of Haule et al. (see Ref. [12]), which is exploited in our numerical implementation.
Given a set of P orbitals, we introduce an orthonormal subset of "uncorrelated" orbitals Q spanning the orthogonal complement to the corresponding P linear space so that P and Q are a complete basis for the physical system considered.
For later convenience, we expand the field operator as where ξ ki;π ðrÞ represents the πth correlated orbital within the i local space (in the Bloch representation), χ σ 0 ðσÞ represents the eigenstate of the third component of the spin with eigenvalue σ 0 , and ΞQ is the component of the field operator that corresponds to the Q orbitals.Equation ( 49) provides a prescription to express any one-body operator A in second quantization as where Let us define V P i as the ith single-particle subspace and Π i as the corresponding orthogonal projector.For later convenience, we also define where GðiωÞ is the Green's function of the system within the whole single-particle space, and ΣðiωÞ is the corresponding self-energy, which is block diagonal in i, within both DMFT and the GA.
The purpose of this section is to derive a functional formulation of the LDA þ GA method [23] with the same mathematical structure of LDA þ DMFT [7].
In Ref. [7], it was shown that the solution of the LDA þ DMFT method can be formulated as the saddle point of the following functional [64]: where μ is the chemical potential, ρðrÞ is the electron density, Δ is the Laplacian, J ðrÞ is the corresponding constraining field, and is the corresponding operator.The functional Φ L i is the i-impurity Luttinger functional associated with a local interaction operator Ĥint i , and Φ dc is an appropriate double-counting correction.In particular, in this work, we assume that Ĥint i is the general (rotationally invariant) Slater-Condon parametrization of the on-site interaction, which is identified by the (atom-dependent) interactionstrength parameters U i and Hund's coupling constant J i , and we employ the following standard form for the doublecounting functional [51]: where is the local population of the i correlated electrons.The total number of electrons N (correlated and uncorrelated) is predetermined by the charge-neutrality condition of the system.
For later convenience, in the rest of this subsection, we assume a linear double-counting functional, represented as where V dc i is a given real number.The generalization to the problem of nonlinear double-counting functionals-such as Eq. ( 57)-will be obtained in Sec.IV D by reducing it to the simpler case of linear double counting.
Note that, under the assumption [Eq.( 59)], the LDA þ DMFT functional can be written as where is the DMFT approximation to the BK functional for the Kohn-Sham-Hubbard (KSH) Hamiltonian where Nloc i is the number operator for the correlated electrons at the site i.

The LDA þ GA functional
In the previous section, we have shown that, in the case of linear double counting [see Eq. ( 59)], the LDA þ DMFT functional can be rewritten as in Eq. ( 60), where Ω KSH V dc ;N is the DMFT functional [Eq.(61)] of the Hubbard model ĤKSH V dc ½J [see Eq. ( 62)].This point of view suggests a natural method to derive the LDA þ GA functional.In fact, our derivation of the LDA þ GA functional consists in replacing in Eq. ( 60) the DMFT functional Ω KSH V dc ;N of the Hubbard model ĤKSH V dc ½J [see Eq. ( 62)] with the corresponding GA functional.Note that the GA functional for a generic Hubbard model was already derived in Sec.II [see Eq. ( 30)].
Specializing Eq. ( 30) to Eq. (62) gives Note that the operators Δ and Ĵ , which appear in the Hubbard Hamiltonian [Eq.( 62)], have been split into their local and nonlocal components consistently with the definitions ( 51) and ( 52).This is because the local and nonlocal operators have to be treated differently within the GA [see, e.g., point (2) in Sec.III A].More precisely, the i local parts of Δ and Ĵ have to be treated together with the interaction Ĥint i , while the nonlocal parts are accounted for within the first term of Eq. ( 63) [see Eq. ( 2) and text below].Note also that, since the Q states are uncorrelated, R is (by construction) a block matrix acting as the identity within the space Q and with blocks R i within the corresponding i correlated spaces.The matrices λ and η are instead zero within the Q space and with blocks λ i and η i , respectively, within the corresponding i correlated spaces.
In summary, the LDA þ GA functional for linear double counting is given by where Ω KSH V dc ;N is given by Eq. ( 63).As in LDA þ DMFT, the total number of electrons N (correlated and uncorrelated) is predetermined by the charge-neutrality condition of the system.

C. Charge self-consistency and KSH Hamiltonian
In this section, we discuss the general structure of our implementation of the charge self-consistent LDA þ GA method in the case of linear double counting.
The stationarity condition of Eq. ( 64) with respect to ρðrÞ is while the stationarity condition with respect to J ðrÞ is where the GA expectation value of the local part of the density operator is not computed from the quasiparticle Green's function but is computed using jΦ i i according to Eqs. ( 10) and (23).Finally, the stationarity condition with respect to the variables i , n 0 i amounts to solving within the GA the KSH Hamiltonian [Eq.( 62)] following the procedure described in Sec.II B.
The generalization of the above numerical procedure to a more general class of double-counting functionals is given in the next section.

D. The LDA þ GA for general double counting
It is useful to generalize the procedure described above to the case when the double-counting functional is a generic nonlinear function of the on-site occupations of the correlated electrons N loc i , as in Eq. (57).In this case, Eq. ( 55) reduces to Ω N ½ρðrÞ; J ðrÞ; G loc ðiωÞ; Σ loc ðiωÞ; μ For later convenience, we promote the average local occupations to independent variables by adding to Eq. ( 67) the following Lagrange-Legendre term: This step enables us to rewrite Eq. ( 67) as follows: Ω N ½N loc ; V dc ; ρðrÞ; J ðrÞ; G loc ðiωÞ; Σ loc ðiωÞ; μ where E dc i is now a function of the new variable N loc i and Ω V dc ;N is the LDA þ DMFT functional valid for the special case of linear double counting [see Eq. ( 60)].Consequently, using Eq. ( 64), we obtain that the LDA þ GA functional is represented as where Ω V dc ;N is the LDA þ GA functional previously derived for the case of linear double counting [see Eq. ( 64)].Equation ( 71) enables us to reduce the LDA þ GA problem for generic double counting to the problem solved in Sec.IV C. In fact, the saddle-point condition with respect to the variables ρðrÞ, J ðrÞ, jΦ i i, and V dc i -of the so-modified Eq. ( 63) can be solved numerically following the procedure of Sec.IV C (see Fig. 1).
The stationarity conditions with respect to N loc i and V dc i give the additional equations which determine the self-consistent V dc i and can be solved numerically, e.g., as shown in Fig. 2.
Note that the structure of the algorithm discussed above and represented in Figs. 1 and 2 In the previous subsection we have shown that the LDA þ GA and the LDA þ DMFT methods require, in order to determine the double-counting potentials V dc i and the charge density ρðrÞ, solving iteratively the correlated Hubbard Hamiltonian [Eq.(62)].
We observe that, remarkably, the calculation of V dc i and ρðrÞ, as well as the total energy, does not require us to compute the spectral properties of the system but only the ground state.
Since the GA ground-state properties are generally in very good agreement with DMFT for strongly correlated metals [23,27,35,65,66]-even though the GA is much less computationally demanding-this observation opens up the possibility to use the GA for structural relaxation and to determine V dc i and ρðrÞ, and to perform a single DMFT iteration afterwards, in order to also have access to the spectral properties of the system of interest.

F. Summary of the main results of Sec. IV
We have derived a functional formulation of the LDA þ GA method [23] with the same mathematical structure of LDA þ DMFT [7] [see Eq. ( 71)].This parallelism has enabled us to use the same LAPW interface between DMFT/GA and the LDA code, and it suggests possible synergistic combinations between the two methods.
We have derived a very stable and numerically efficient implementation of the LDA þ GA method, whose structure is as follows.(1) The double-counting potentials V dc i are determined by the outer loop represented in Fig. 2. (2) Each iteration of the outer loop requires us to calculate the LDA þ GA solution at fixed V dc i -i.e., the saddle point of the functional [Eq.( 63)]-which is computed numerically using the charge self-consistent procedure represented in Fig. 1.Each iteration of the charge self-consistency loop consists in solving the Kohn-Sham-Hubbard Hamiltonian determined by the input electron-density according to Eq. ( 62) and computing the corresponding output electron density according to Eq. ( 66) until convergence.(3) The Kohn-Sham-Hubbard Hamiltonian is solved using the procedure given in Sec.II B.

V. THE PRESSURE-VOLUME PHASE DIAGRAM OF PRASEODYMIUM
In this section, we apply our LDA þ GA implementation to the elemental praseodymium.As in Refs.[35,36], we employ the "standard" prescription for the double-counting functional and the general Slater-Condon parametrization of the on-site interaction [51], assuming that the Hund's coupling constant is J ¼ 0.7 eV and that the value of the interaction strength is U ¼ 6 eV, which is consistent with previous constrained LDA calculations [50].However, since the value of U is generally difficult to predict exactly, here we also perform calculations for U ¼ 5 eV and U ¼ 7 eV.
The lattice structure of the elemental praseodymium is dhcp at ambient conditions, and it undergoes the following sequence of transformations under pressure: dhcp → fcc → distorted − fcc → α-U [50,67].While no appreciable volume collapse is associated with the transitions between the three lower-pressure phases-which are all characterized by relatively high symmetry and/or good packing ratios [68]-the distorted − fcc → α-U transition is accompanied by a sizable volume collapse (about 10% at room temperature).
It is widely believed that the low-symmetry α-U lattice structure in the elemental praseodymium is stabilized at high pressures by the delocalized f electrons, in accordance with a general argument based on the Peierls theorem FIG. 2. Schematic flow chart of the LDA þ GA outer loop that determines the double-counting potentials V dc i .The step concerning the solution of the problem at fixed V dc i is implemented as in Fig. 1 (see Sec. IV C).
(see Refs. [69,70]).The deformation from a high-symmetry structure to a low-symmetry structure can lower the band's energy by opening a Peierls gap between a "bonding" band (below the Fermi level) and an "antibonding" band (above the Fermi level).This effect competes with the electrostatic Madelung interaction, which favors the high-symmetry lattice structures [71] such as the fcc.Nevertheless, there are a couple of important aspects of the physics underlying the volume-collapse transition of Pr that are still not fully understood and that require further investigation.( 1) What is the role of the electron correlations for the volume-collapse transition?(2) Would Pr display a volume-collapse transition even without taking into account the change of structure?
In order to further investigate the physics underlying the volume-collapse transition of Pr, in this work, we study it theoretically in the fcc and α-U lattice structures (see Fig. 3).
The theoretical LDA and LDA þ GA phase diagrams are shown in the lower panels of Fig. 4, in comparison with the experimental data.The pressure was obtained as P ¼ −dE=dV from the corresponding total-energy curves, which are shown in the upper panels.Note that the calculations at U ¼ 6 eV are shown in the main panels, while the calculations at U ¼ 5 eV and U ¼ 7 eV are shown in the insets.The LDA þ DMFT calculations of fcc Pr of Ref. [72], which were performed at T ¼ 632 K, are also reported.Note that in these LDA þ DMFT calculations the charge self-consistency was not carried out, and part of the Slater integrals and the non-density-density terms in the local f-electron interaction were neglected.The theoretical and experimental values of equilibrium volume, bulk modulus, and critical pressure are reported in Table I.
The agreement between the LDA þ GA theoretical results and the experimental data is very good.In particular, U ¼ 7 eV gives an overall better agreement with the experimental phase diagram.Remarkably, the GA correction to the LDA total energy is very important for Pr, especially at low pressures.In particular, the LDA equilibrium point is about V LDA eq ≃ 21 Å 3 =atom, while the experimental value is approximately V exp eq ≃34.5Å 3 =atom, which is better reproduced by the LDA þ GA calculations, which give V LDAþGA eq ≃ 32 Å 3 =atom.Furthermore, while FIG. 4. Theoretical total energy as a function of the volume (upper panels) at zero temperature for the fcc and the α-U phases.The corresponding pressure-volume phase diagrams (lower panels) are shown in comparison with the experimental data of Ref. [50] (colored markers) and Ref. [67] (gray markers), which refer to measurements at room temperature.The LDA calculations are reported in the left panels, and the LDA þ GA calculations are reported in the right panels.The LDA þ DMFT calculations of fcc Pr at T ¼ 632 K of Ref. [72] are also shown (see the dotted red lines in the right panels). the LDA predicts that the α-U structure becomes less stable than the fcc phase only at negative pressures, the LDA þ GA method predicts correctly that the volume collapse occurs at positive pressures (see the commontangent construction in the upper-right panel of Fig. 4).Our LDA þ GA calculations are also consistent with the LDA þ DMFT calculations of fcc Pr of Ref. [72].
Our theoretical results indicate that in Pr, the correlation effects energetically favor the fcc lattice structure with respect to the α-U, stabilizing it at the equilibrium point and for a wide range of positive pressures.This fact is also clearly illustrated in the inset of the upper-left panel of Fig. 4, which represents the energy difference between the fcc and α-U structures as a function of the volume, both in LDA and in LDA þ GA.
We point out that, contrarily to Ce [35], the fcc phase of Pr would not display any isostructural volume-collapse transition by applying pressure.In fact, the second derivative of the fcc energy-volume curve is positive within the entire range of volumes and for all of the values of U considered.
Let us examine how the correlation effects taken into account by the GA correction influence the on-site f occupation probabilities W f and the f quasiparticle renormalization weights, which are determined as Z ¼ R † R according to Eq. ( 8).In the upper panels of Fig. 5, the on-site f occupation probabilities are illustrated for the fcc and for the α-U phases as a function of the volume.While the LDA probability distribution is very "broad" at all pressures, the LDA þ GA probability distribution is relatively narrow, especially at large volumes, where we find that the majority of the f electrons lie within the f 2 space, in agreement with recent experiments [69].The averaged f quasiparticle renormalization weights are shown in the upper panels of Fig. 6.Note that, because of the spin-orbit coupling, the f quasiparticle weights are split into two groups with total angular momentum J ¼ 5=2 and 7=2, respectively.While at small volumes both of the Z's decrease as a function of the volume, at larger volumes they develop a qualitatively different behavior: Z 5=2 becomes significantly smaller than 1-indicating that the system is very correlated in this regime-and Z 7=2 increases.This behavior is a consequence of the spin-orbit FIG. 6. Evolution as a function of the volume of the averaged LDAþGA (U ¼ 6 eV) quasiparticle renormalization weights Z (upper panels) of the averaged 5=2 and 7=2 f electrons and corresponding orbital occupations (lower panels) for the fcc phase (left panels) and for the α-U phase (right panels) of Pr.For the fcc phase, the LDAþDMFT results of Ref. [72] are also reported (dots).FIG. 5. Theoretical on-site f occupation probabilities W f for the fcc phase (left panels) and for the α-U phase (right panels).The LDA results are reported in the upper panels, and the LDA þ GA results at U ¼ 6 eV are reported in the lower panels.For the fcc phase, the LDA þ DMFT results of Ref. [73] are also reported (dots).effect, which occurs also in Ce [35,36].In particular, Z 7=2 grows because the 7=2 electrons disappear at larger volumes, as indicated in the upper panels of Fig. 6.Note that the difference for W f and the Z's is very small between the two lattice structures.
As shown in the lower panels of Fig. 4, the GA correction to the pressure-volume phase diagram becomes more substantial at large volumes, where the quasiparticle renormalization weights are considerably smaller than 1.This is not surprising, as decreasing the distance between the atoms increases the bandwidth of the f electrons, which reduces the relative importance of the GA correction to the LDA functional-i.e., of the interaction and the doublecounting terms [see Eq. ( 63)].Nevertheless, it is important to observe that the above-mentioned correction to the pressure-volume phase diagram is essentially identical for the two lattice structures considered.Note that this consideration is consistent with the fact that, according to our calculations, the strength of the electron correlations is very similar for the fcc and the α-U lattices.
Let us now analyze the role of the electron correlations in the determination of the more stable structure as a function of the volume.As we have anticipated, the Peierls mechanism, which stabilizes the α-U phase, relates to the itinerant character of the f electrons, and it is consequently less effective at large volumes, where the f bandwidth is smaller.Note that this effect is qualitatively well captured already by the LDA.In fact, at very large volumes (negative pressures), the LDA total energy of the fcc lattice structure becomes lower with respect to the α-U (see the upper-left panel of Fig. 4).As shown in the inset of the upper-right panel of Fig. 4, the effect of the GA correction on the total energy difference between the two lattice structures is of the same order of magnitude for all volumes.Consequently, we attribute the consequent improved quantitative agreement with the experiments for the transition volume to the overall more realistic evaluation of the total energy in LDA þ GA with respect to LDA.In other words, we argue that the behavior of the energy difference between the fcc and the α-U structures is essentially already qualitatively captured by the LDA, and it is not directly related to the f-electron localization-which is substantial only at large volumes, as indicated by the f quasiparticle weights shown in the lower panels of Fig. 6.
In conclusion, we have observed that the behavior of the quasiparticle weights and of the f configuration probabilities as a function of the volume is essentially equal for the fcc and the α-U lattice structures.In particular, at small volumes we find that Z ≃ 1, while at large volumes we find that Z ≪ 1. Consistently with the Peierls mechanism, while at small volumes the α-U structure has a lower energy, at large volumes the fcc lattice configuration becomes more stable.Since both the energy difference between the phases and the quasiparticle weights are controlled by the volume, it is not surprising that the volume-collapse transition of Pr is accompanied by an abrupt delocalization of the f electrons, and it is not even surprising that a correlation between these two phenomena is found experimentally in several other rare-earth materials [50].On the other hand, based on our calculations, it does not seem appropriate to regard the above-mentioned correlation as a general causeeffect relation.In fact, Pr would not display the transition without taking into account the change of structure (at least at low temperatures).Note also that other f systems, such as americium [74], display volume-collapse structure transitions maintaining essentially a constant f valence, indicating that f localization is not a crucial prerequisite for the volume-collapse transitions in f systems.

VI. PHASE DIAGRAM OF PLUTONIUM
Plutonium is the most exotic and mysterious element in the periodic table.Its stable structure at ambient conditions is α-Pu, which has a low-symmetry monoclinic structure with 16 atoms within the unit cell grouped in eight nonequivalent types.At higher temperatures (see Fig. 7), Pu can assume the following distinct lattice structures: β (monoclinic, with 34 atoms within the unit cell grouped in seven inequivalent types), γ (orthorhombic), δ (fcc), δ 0 (bct), and ϵ (bcc).One of the most intriguing properties of Pu is that these temperature-induced structure transitions are accompanied by significant changes of density.In particular, the equilibrium volumes of the δ and δ 0 phases are very large with respect to the other allotropes.Another interesting puzzle is that δand α-Pu have negative thermal-expansion coefficients within their respective range of stability, unlike the vast majority of materials.These facts have stimulated extensive theoretical and experimental studies.However, a convincing explanation of the metallurgic properties of Pu based on fundamental principles is still lacking, and none of the previous theories has been able to describe simultaneously the energetics and the f electronic structure of all of the phases of Pu on the same footing.
Previous state-of-the-art DFT calculations [75][76][77][78] were able to reproduce reasonably well the equilibrium volumes of the phases of Pu, but in order to describe all of them on the same footing, it was necessary to introduce artificial [79] spin and/or orbital polarizations-thus compromising the description of the electronic structure.In fact, without spin and orbital polarization, these techniques predict that the equilibrium volumes of all of the phases are essentially identical, in contrast with the experiments (see Fig. 7).Calculations within the framework of DFT in combination with dynamical mean field theory (DFT þ DMFT) have explained several aspects of the electronic structure of Pu (see, e.g., Refs.[80][81][82][83][84]).Nevertheless, the computational complexity of this approach made it impossible to calculate the pressure-volume phase diagram of all of the phases of Pu.
In this section, we provide a bird's eye view of Pu by studying all of its crystalline phases at zero temperature using our implementation of the LDA þ GA method, whose description of the ground-state properties is generally in very good agreement with LDA þ DMFT but is considerably less computationally demanding.In particular, we employ the general Slater-Condon parametrization of the on-site interaction with parameters U ¼ 4.5 eV and J ¼ 0.36 eV-as we find that LDA þ GA calculations performed using these values give a better overall agreement with the experiments with respect to U ¼ 4.5 eV and J ¼ 0.51 eV (which are the values previously assumed in Ref. [84]).Calculations of α-Pu with different values of U and J are shown in Ref. [85].Very interestingly, our study indicates that the electron correlations are only weakly dependent on the lattice structure, while the most important element originating the differentiation between the equilibrium densities of the phases of Pu is the competition between the Peierls effect and the Madelung interaction [70,71].
In the upper panels of Fig. 8, we show the LDA (left) and LDA þ GA (right) evolutions of the total energies E as a FIG. 8. Theoretical total energies for the crystalline phases of Pu as a function of the volume (upper panels) and corresponding pressure-volume curves (lower panels) in comparison with the experimental data of α-Pu from Refs.[88] (black circles), [89] (blue squares), and [86] (red diamonds).Our results are shown both in LDA (left panels) and in LDA þ GA (right panels).The right insets are zooms of the curves in the corresponding panels.In the upper-left inset, we show the correlation energies, while the corresponding contributions to the pressure are shown in the lower-left inset.The vertical lines indicate the minima of the energy curves.The horizontal lines of the legend indicate the estimated zero-temperature equilibrium volumes, which are assumed to lie between the zero-temperature values extrapolated by linear interpolation and the experimental values at the temperatures in which the allotropes are stable (see Fig. 7).function of the volume V for all of the crystalline phases of Pu.In the lower panels, we show the corresponding evolutions of the pressure P ¼ −dE=dV in comparison with the experimental data of α-Pu.In the left insets, we show the correlation energies-here defined as the differences between the LDA þ GA and LDA total energies-and the respective contributions to the pressure.Note that in our calculations, we have not performed structure relaxation, but we have assumed a uniform rescaling of the experimental lattice parameters (see Refs. [86,87]).
In Table II, the theoretical zero-temperature equilibrium volumes are shown in comparison with the zerotemperature experimental volumes, which we assume to be in between the thermal-equilibrium volumes and the zero-temperature values extrapolated by linear interpolation in Fig. 7.The bulk modulus and energies (referred to as the ground-state energy of α-Pu) are shown in comparison with the experimental data of Refs.[78,90].Remarkably, while the theoretical equilibrium volumes of all phases of Pu are very similar in LDA, they are very different in LDA þ GA, and in good quantitative agreement with the zero-temperature experimental values.Furthermore, while LDA predicts very large equilibrium energy differences between the phases of Pu, these differences are very small in LDA þ GA, in agreement with the experiments.Note also that the LDA þ GA ground-state energies increase monotonically from each phase to the next-highertemperature phase, consistently with the experiments (see Fig. 7 and Table II).The only exception is β-Pu, whose theoretical equilibrium energy is larger than γ-, δ-, and δ 0 -Pu.
In order to understand how the electron correlations so drastically affect the energetics of Pu, it is enlightening to look at the behavior of the correlation energies (see the left insets in Fig. 8).In fact, the evolution of the correlation energies as a function of the volume is essentially structureless and identical for all of the phases (except for a uniform structure-dependent energy shift whose main effect is to slightly increase the energy of α-Pu with respect to the other phases).As a result of this correction, the LDA total energies are transformed as indicated by the gray circles in the upper panels of Fig. 8.The relative behavior of the LDA þ GA zero-temperature energies of the allotropes of Pu is clearly inherited by the LDA energy-volume curves in the region highlighted in the upper-left panel of Fig. 8, which transform into the region highlighted in the upper-right panel of Fig. 8 when the correlation energies are taken into account.The same considerations apply to the evolutions of the pressure, as indicated by the gray circles in the lower panels of Fig. 8.
The above observations explain from a simple perspective how the electron correlations determine the unusual energetics of Pu.In fact, the energy crossings between the LDA energy curves of the high-symmetry structures (δ-and α-Pu) and the other phases can be simply understood in terms of the competition between the Peierls effect and the Madelung interaction-which is known to energetically favor the low-symmetry structures at small volumes and the high-symmetry structures at large volumes [69][70][71].The most important effect of the correlation energies is to shift the position of the equilibrium volumes to larger values, near where the above-mentioned energy crossings take place and the LDA energy differences are relatively small.Interestingly, a similar interplay between correlation effects and bands structure is displayed in Pr (see Sec. V) and might emerge in even greater generality.
In the upper panels of Fig. 9, we show the occupations of the f electrons.The total number of f electrons in δ-Pu is n f ≃ 5.2, which is consistent with previous LDA þ DMFT calculations [81].Here, we find that n f ≃ 5.2 also for γ-, δ 0 -, and ϵ-Pu.In the monoclinic structures, n f is different for the inequivalent atoms within the unit cell; it runs between 5.21 and 5.32 in α-Pu, while it runs between 5.17 and 5.21 in β-Pu.In the middle panels of Fig. 9 we show the averaged orbital populations with total angular momentum J ¼ 7=2 and J ¼ 5=2.Note that for all of the phases of Pu, the number of 7=2 f electrons decreases as a function of the volume, while the number of 5=2 f electrons increases.This behavior simply indicates that the spin-orbit effect is more effective at larger volumes, as expected.Finally, in the lower panels, we show the behavior of the branching ratio B, which is a measure of the strength of the spin-orbit coupling interaction in the f shell and is calculated from the orbital populations using the following equation: TABLE II.Zero-temperature theoretical equilibrium volumes, bulk modulus, and total energies of the crystalline phases of Pu in comparison with the experiments [78,90].Consistently with Fig. 8, the zero-temperature equilibrium volumes are assumed to lie between the zero-temperature values extrapolated by linear interpolation and the experimental values at the temperatures in which the allotropes are stable (see Fig. 7).

Pu
(see Refs. [91,92]).Consistently with the behavior of the orbital populations, B increases as a function of the volume.Note that the behavior of n f and B is very similar for all of the phases of Pu.
The α-Pu theoretical value of n f at equilibrium is in good agreement with the values extrapolated from the x-ray absorption near-edge structure (XANES) measurements of Ref. [93].On the other hand, while our calculations indicate that n f is slightly smaller in δ-Pu than in α-Pu, according to the extrapolations of Ref. [93], the 1.9%-Ga δ-Pu alloy has a larger n f with respect to α-Pu.Also, the theoretical values of B are in good agreement with values extrapolated in Refs.[94,95] from electron energy-loss spectroscopy (EELS) and x-ray absorption spectroscopy (XAS) [91,[96][97][98].
Let us study the behavior of the many-body reduced density matrix ρf of the f electrons, which is obtained from the full many-body density matrix of the system by tracing out all of the degrees of freedom with the exception of the f local many-body configurations of one of the Pu atoms.We define where k is an arbitrary constant that we determine so that the lowest eigenvalue of F is zero by definition.Within this definition, ρf ∝ e − F; i.e., F represents an effective local Hamiltonian of the f electrons that depends on the volume and that is renormalized with respect to the atomic f Hamiltonian because of the entanglement with the rest of the system (see Ref. [99]). Figure 10 shows the eigenvalues P n of ρf as a function of the eigenvalues f n of F for all of the allotropes of Pu (that are computed at their respective theoretical zerotemperature equilibrium volumes).Consistently with Ref. [81], we find that for δ-Pu there are two dominant groups of multiplets: one with N ¼ 5 and J ¼ 5=2 (that is 6 times degenerate) and one with N ¼ 6 and J ¼ 0 (that is nondegenerate).Interestingly, our results show that this conclusion also applies to all of the other phases.The f probability distribution of δand δ 0 -Pu is slightly less broad FIG.10.Configuration probabilities of the eigenstates of the reduced density matrix ρf ≡ e − F=Tr½e − F of the f electrons as a function of the eigenvalues f n of F. Each configuration probability is weighted by the degeneracy d n ¼ 2J n þ 1 of the respective eigenvalue f n , where J n is the total angular momentum.FIG. 9. Upper panels: Evolution as a function of the volume of the averaged orbital populations of the 5=2 and 7=2 f electrons.Middle panels: Total orbital occupations in comparison with the values extrapolated in Ref. [93] from XANES measurements at ambient conditions of α-Pu (black stars) and the 1.9%-Ga δ-Pu alloy.Lower panels: Theoretical branching ratios in comparison with the values extrapolated in Refs.[94,95] from XAS (black cross) and EELS (blue crosses) experiments of α-Pu and the 0.6%-Ga δ-Pu alloy.The colors in the first and second panels from the left correspond to the inequivalent atoms of α-Pu and β-Pu.The vertical dotted lines indicate the LDA þ GA equilibrium volumes for the respective phases.
with respect to the other phases.This is to be expected, as δand δ 0 -Pu are stable at larger volumes, such that the local f degrees of freedom are less entangled with the rest of the system with respect to the other phases.The f probability distributions of αand β-Pu are considerably different for inequivalent atoms, as they depend on the number and relative distances of the nearest-neighbor atomic positions, consistently with Ref. [84] (see also Ref. [85]).Note that in β-Pu the atom dependency of the f probability distribution is less pronounced than in α-Pu.
We point out that n f ≃ 5.2 reveals that the f electrons of Pu are in a pronounced mixed-valence state [100].Indeed, the probability of the N ¼ 6, J ¼ 0 multiplet is very large, as indicated by the fact that it has the lowest F eigenvalue.The reason why n f is closer to 5 than to 6 is that the N ¼ 6, J ¼ 0 multiplet is nondegenerate, while the N ¼ 5, J ¼ 5=2 F eigenvalue is 6 times degenerate-so its contribution to n f is weighted by a factor 6 ¼ 2 × 5=2 þ 1.The observation that the f electrons have a significant mixed-valence character indicates that the local f degrees of freedom are highly entangled with the rest of the system [99].This observation is consistent with the fact that the Pauli susceptibility of the δ-Pu Ga alloy is Pauli-like at low temperatures [79].Furthermore, it is consistent with the statement of Ref. [90] that Pu is an ordinary quasiharmonic crystal in all of its crystalline phases; i.e., already at T ≳ 200 K, the electronic entropy is very small with respect to the quasiharmonic contributions.
In conclusion, we have calculated from first principles the zero-temperature energetics of Pu, finding good agreement with the experiments.Our analysis has clarified how the electron correlations determine the unusual energetics of Pu, including the fact that the different allotropes have very large equilibrium-volume differences while they are very close in energy.Remarkably, in our calculations, we did not introduce any artificial spin and/or orbital polarizations [79], while this was necessary in previous state-ofthe-art DFT calculations [75][76][77][78].This advancement has also enabled us to describe the f electronic structure of Pu on the same footing.Our calculations indicate thatsimilarly to Pr-the ground-state f electronic structure is similar for all of the phases of Pu and that the f-electron atomic probabilities display a significant mixed-valence character.Our zero-temperature calculations of Pu also constitute an important step toward the theoretical understanding of its peculiar temperature-dependent properties, e.g., the negative thermal expansion of δand α-Pu.In fact, above room temperature, the contributions to the free energy of the nonadiabatic effects and of the thermal excitations of the electrons from their ground state are expected to be negligible in Pu [90].Consequently, the total energy could be used to investigate this problem either within direct Monte Carlo simulations or using molecular dynamics [101,102] with atomistic potentials extrapolated from this work.

VII. CONCLUSIONS
We have developed an exceptionally efficient algorithm to implement the GA.Furthermore, we have derived a functional formulation of LDA þ GA that has the same mathematical structure of LDA þ DMFT [7].This insight has enabled us to pattern the LAPW interface [44] between LDA and GA after the LDA þ DMFT work of Ref. [12].Using our LDA þ GA code, we have performed firstprinciples calculations of Pr and Pu under pressure, which are prototypical systems with partially localized f electrons.
Our calculations of Pr indicate that its volume-collapse transition is not driven only by the concomitant delocalization of the f electrons.In fact, contrarily to Ce [35], Pr would not display any volume-collapse transition without taking into account the change of lattice structure (at least at low temperatures).This suggests that there is no reason to exclude the possibility that, in other f materials, a volumecollapse transition may occur without any concomitant substantial f delocalization-as indicated, for instance, by recent experiments on the elemental Am [74].Note that the understanding of the connection between f delocalization and volume-collapse transitions in f systems is one of the most important puzzles in condensed matter theory.
Our calculations of Pu constitute the first theoretical description of all of the crystalline phases of this material giving good agreement with all of the experiments on the same footing-including both the thermodynamical properties and the f electronic structure.A particularly important conclusion of our study is that the most important effect originating the differentiation between the equilibrium densities of the phases of Pu is the competition between the Peierls effect and the Madelung interaction and not the dependence of the electron correlations on the lattice structure, which is a negligible effect.Note that the explanation of this phenomenon is of great interest both physically and from a metallurgic standpoint.
From a technical point of view, our calculations clearly demonstrate the exceptional capabilities of the computational scheme presented in this work.Indeed, our method enables us to rapidly perform accurate first-principles calculations of strongly correlated materials even for systems so complex that other state-of-the-art methods are too time-consuming to be practically applicable.
Department of Energy by Iowa State University under Contract No. DE-AC02-07CH11358.

APPENDIX A: SUMMARY OF THE STANDARD FORMULATION OF THE GUTZWILLER VARIATIONAL METHOD
Let us write the Hubbard Hamiltonian in an infinitecoordination lattice [see Eq. ( 1)] in the form where here, in order to derive our formalism, we conveniently assume that jΓ; Rii are local Fock states, From now on, we will name the above c basis set the "original" basis.The GA consists in variationally determining a projected wave function represented as where jΨ 0 i is a Slater determinant and PRi is a general operator acting on the local configurations at the site ðR; iÞ, which we represent in the original basis as follows: Here, we assume that jΨ 0 i is an eigenstate of the number operator and that the local operators PRi satisfy the commutation rule PRi ; which implies that jΨ G i is also an eigenstate of the number operator.Note that these assumptions are no longer valid when the method is generalized to study superconductivity (see, e.g., Refs.[22,25]), but this subject will not be addressed in the present work.
In order to analytically evaluate the expectation value of Ĥ with respect to jΨ G i, the following two additional approximations are done.
(1) The manifold of variational wave functions is further restricted by the following conditions [22]: which are commonly named "Gutzwiller constraints." (2) The so-called GA is assumed, which is an approximation scheme that, like DMFT [9], becomes exact in the limit of infinite coordination lattices (see, e.g., Ref. [22]).
The derivation of the Lagrange formulation of the GA [Eq.( 3)]-that we employed in Sec.II as a starting pointis summarized in this appendix.The material presented in this appendix makes large use of ideas developed in previous works by several authors [22][23][24][25]27,103].

Reformulation of the Gutzwiller problem
In this section, we briefly summarize the reformulation of the Gutzwiller problem derived in Refs.[24,25].

a. The mixed-basis representation
Let us consider the matrix Since ρ 0 i are Hermitian, there always exists a unitary transformation U i that diagonalizes it, i.e., such that The so-obtained ladder operators f Rib are named naturalbasis [22] operators.
Note that the connection between the natural basis f and the original basis c depends only on jΨ 0 i.At given jΨ 0 i, the coefficients Λ that determine the Gutzwiller projector [see Eq. (A4)] are free variational parameters.
Instead of expressing the Gutzwiller projector in terms of the original basis as in Eq. (A4), it is convenient to adopt the following mixed original-natural [24] For later convenience, we adopt the convention that the order of the jΓ; Rii and the jn; Rii states is the same.For instance, if the second Γ vector in Eq. (A12) is ĉ † 1↑ ĉ2↓ j0i, then the second n vector is f † As we will see, the mixed-basis parametrization of the Gutzwiller projector enables us to gauge away from the formalism the unitary matrix U that relates the original basis and the natural basis [24] [see Eq. (A10)], which is a great simplification.

b. Gutzwiller expectation values
In an infinite-coordination lattice, the expectation value of any observable can be computed analytically.
Let us define the uncorrelated occupation-probability matrices P 0 i with elements [22] ½P 0 i nn 0 ≡ hΨ 0 jjn 0 ; Riihn; RijjΨ 0 i We also introduce the matrix representations of the operators fRib and ĉRiβ , Note that, since we have assumed that the order of the jΓ; Rii and the jn; Rii states is the same, the matrix elements of F ib can be equivalently computed as With the above definitions, it can be readily verified that the expectation value of any local observable can be calculated as where and that the Gutzwiller constraints can be written as The average of the intersite density matrix reduces to [24,26] In other words, the intersite single-particle density matrix averaged on jΨ G i is computed by averaging over jΨ 0 i a renormalized density matrix with natural fermionic operators, replacing the physical ones according to the rule Within the definitions given in this section, it can be shown that the renormalization matrices R i can be expressed as Note that, thanks to the mixed-basis representation of the Gutzwiller projector, the unitary transformation that relates the natural-basis operators f to the original ones c need not be known explicitly, which is a great simplification.

c. The ϕ matrix and the GA total energy
The formalism can be further simplified by defining the following matrix [24]: Note that Eq. (A5) translates into Eq.( 6) for ϕ i .Within this definition, the expectation value of any local observable [see Eq. (A17)] reduces to the Gutzwiller constraints [see Eqs.(A19) and (A20)] simplify as follows: and Eq.(A23) for the renormalization factors reduces to In conclusion, the variational energy [Eq.(A6)] is given by [24] and it has to be minimized by fulfilling Eqs.(A26)-(A28).Following Ref. [27], we take into account the constraints by applying the theorem of the Lagrange multipliers as follows: Note that n 0 i and R i have been conveniently promoted to independent variables [104].Within this formulation, the numerical problem arising from the stationarity condition of L is particularly easy to solve numerically (see Sec. II B).
It is convenient to rewrite Eq. (A30) as follows: where is the GA quasiparticle Hamiltonian of the system [105].

d. GA Lagrange functional in the canonical ensemble
If the Hubbard Hamiltonian [see Eq. (A1)] has to be solved in the canonical ensemble, i.e., at a fixed number of particles per site N, the functional [Eq.(A31)] becomes where Note that the constraint on the number of particles has been imposed on jΨ 0 i instead of jΨ G i.This is licit because we have assumed that the Gutzwiller projector PG commutes with the number operator N. where we used N n 0 ¼ N n þ 1, so e −iðπ=2ÞN n 0 ðN n 0 −1Þ e iðπ=2ÞN n ðN n −1Þ ¼ e −iπN n ¼ ð−1Þ N n : ðB18Þ It can be readily verified that where we introduced the matrix elements of F † α , which is the representation of the ladder operators in their own Fock basis.The phase ð−1Þ N n appearing in the first step of Eq. (B19) is due to the fact that the c and f ladder operators anticommute, so a minus sign is generated whenever they are transposed.Note that in the last step we used that F † a is real (its matrix elements can only assume the values 0, 1, and −1), so doing the Hermitian conjugate amounts to transposing the labels.

APPENDIX C: EXPECTATION VALUES IN TERMS OF THE LOCAL GREEN'S FUNCTIONS
In this section, we derive useful alternative expressions for Eqs. ( 32)- (35).
It is well known that the Fermi function can be expressed as a partial fraction decomposition (PFD), where several alternative choices are possible for z n .For instance, the well-known Matsubara expansion is given in terms of the Matsubara frequencies z n ¼ iω n , which are purely imaginary.However, the convergence of the Matsubara series is very slow.A numerically more convenient alternative to the Matsubara expansion was proposed in Ref. [106], which converges faster than exponentially, and where the poles z n are complex, with finite real and imaginary components.Equation [106] enables us to express the left members of Eqs. ( 32)- (35) in terms of the local quasiparticle Green's functions as follows: We observe that G qp i can be rewritten as PHASE DIAGRAM AND ELECTRONIC STRUCTURE OF … PHYS.REV.X 5, 011008 (2015) 011008-23 where ΣðzÞ, which represents the GA approximation for the self-energy, is a block-matrix Eq. ( 8)] whose blocks are given by and is the GA approximation for the coherent part of the local Green's function.Within the above definitions, it can be readily shown that Eqs.(C3) and (C4) can also be represented as Note that the evaluation of Eqs.(C8) and (C9) only requires knowledge of the local Green's functions G i of the correlated sites [see Eq. (C7)].As we will show, the operation to compute G i at given Σ i relates to the embedding procedure of DMFT [9].We observe that the computational time necessary to calculate Eq. (C7) scales as the square of the number of correlated atoms.On the contrary, Eqs. ( 32)-( 35) require the diagonalization of ϵ k , whose computational time scales as the cube of the number of total atoms (correlated and not).Consequently, Eqs.(C8) and (C9) can be preferable if the unit cell contains many atoms and/or if only a few orbitals are correlated.

APPENDIX D: APPLICATION TO IMPURITY MODELS
In the previous section, we discussed our method to solve the GA equations for a generic Hubbard model.Let us now consider a generic impurity Anderson model (IAM), where Ĥloc is the Hamiltonian of the impurity, which includes both the interaction and the one-body component.The purpose of this section is to solve ĥ within the GA in the grand-canonical ensemble.
We observe that any IAM can be viewed as a Hubbard model, with no translational invariance and where only the impurity site is correlated.In this special case, the GA equations derived in the previous section reduce to where is the hybridization function, is the coherent part of the impurity Green's function, and is the Gutzwiller self-energy of the impurity (that also includes the on-site energies in our notation).Note that the impurity quasiparticle Green's function that corresponds to Eq. (D11) is given by

FIG. 1 .
FIG. 1. Schematic flow chart of the LDA þ GA charge selfconsistent procedure for linear double counting.The solution of ĤKSH is calculated as discussed in Sec.II B.

FIG. 3 .
FIG. 3. Representation of the fcc and α-U crystal structures of Pr.

FIG. 7 .
FIG. 7. Experimental volume-temperature phase diagram of Pu.The dotted lines indicate the zero-temperature equilibrium volumes extrapolated by linear interpolation.

TABLE I .
Theoretical equilibrium volume V eq , bulk modulus K, and critical pressure P c of Pr in comparison with the experiments.The LDA results are also given for α-U Pr, which has a lower energy with respect to fcc Pr (within LDA).
basis form