Density functional theory plus dynamical mean field theory within the framework of linear combination of numerical atomic orbitals: Formulation and benchmarks

The combination of density functional theory with dynamical mean-field theory (DFT+DMFT) has become a powerful first-principles approach to tackle strongly correlated materials in condensed matter physics. The wide use of this approach relies on robust and easy-to-use implementations, and its implementation in various numerical frameworks will increase its applicability on the one hand and help crosscheck the validity of the obtained results on the other. In the work, we develop a formalism within the linear combination of numerical atomic orbital (NAO) basis set framework, which allows for merging NAO-based DFT codes with DMFT quantum impurity solvers. The formalism is implemented by interfacing two NAO-based DFT codes with three DMFT impurity solvers, and its validity is testified by benchmark calculations for a wide range of strongly correlated materials, including 3\textit{d} transition metal compounds, lanthanides, and actinides. Our work not only enables DFT+DMFT calculations using popular and rapidly developing NAO-based DFT code packages, but also facilitates the combination of more advanced beyond-DFT methodologies available in this codes with the DMFT machinery.


I. INTRODUCTION
The understanding of strong correlations among electrons in realistic materials is of great importance for both fundamental science and technological applications. The strong electron correlations trigger a variety of exotic phenomena, which provide a fruitful avenue for designing novel materials [1][2][3]. Kohn-Sham (KS) density functional theory (DFT) [4,5] under localdensity approximation (LDA) and generalized gradient approximation (GGA) has achieved remarkable success in describing a very wide range of materials. However, these approximations are found to be inadequate for correctly describing strongly correlated materials, e.g., 3d transitional metal compounds, lanthanides, actinides, quantitatively or even qualitatively [6][7][8][9]. To overcome this limitation, the so-called DFT+DMFT method that merges DFT and many-body technique based dynamical mean-field theory (DMFT) has been developed [7][8][9][10][11][12][13][14][15]. Various studies have shown the prominent strength of DFT+DMFT in describing electronic structures of strongly correlated materials. DFT+DMFT has thus becoming a promising method of choice in studying realistic strongly correlated materials with certain predictive power. [8,9,[16][17][18][19] Comparing to DFT+U [6,20,21], a method that shares similar spirit as DFT+DMFT and is available in nearly all popular DFT code packages, DFT+DMFT is not only computationally more expensive, but also relatively more difficult to access for normal users. Implementing DFT+DMFT with well-defined local orbital basis and in a way that doesn't requires much expertise from users is, thus, of practical importance for promoting the methodology in studying realistic correlated materials. One of the key issues in the implementation of DFT+DMFT scheme is the definition of local correlated orbitals, in terms of which the DMFT impurity problem is defined and solved. Many-body corrections to the KS Hamiltonian exist only in this subspace. As a result, the choice of local correlated orbitals has a noticeable influence on the obtained DFT+DMFT results at the quantitative level.
In early days, DFT+DMFT was implemented within the linear muffin-tin orbital (LMTO) basis set framework [14-16, 22, 23], where the local LMTOs were chosen as the basis orbitals hosting strongly correlated d or f electrons. Later, different Wannier-type orbitals, such as projective Wannier functions [24] and maximally localized Wannier functions, were used to define the local orbitals [25,26]. As for plane-wave based DFT, Amadon et al. [27,28] used the all-electron atomic partial waves within the projector augmented wave (PAW) framework or pseudoatomic wave functions as the local correlated orbitals. Within the linearized augmented plane-wave (LAPW) method, Aichhorn et al. [29], and Haule et al. and Kim [30] independently developed projector schemes that convert the KS orbital space to local correlated subspace. We note that, in recent years, the underlying embedding idea of DMFT has been used in broader sense arXiv:2205.03620v1 [cond-mat.str-el] 7 May 2022 whereby the embedded cluster is solved in a ab initio way using beyond-DFT approaches such as GW or coupledcluster method, whereas the environment is treated using less advanced approaches [31][32][33].
In recent years, numerical atomic orbitals (NAOs) have been employed as the basis set choice in several firstprinciples software packages [34][35][36][37][38]. Unlike most other basis sets, the linear combination of NAOs represents a versatile framework that can be used in both allelectron and pseudopotential-based DFT calculations. Past experience indicates that NAOs are an efficient basis set choice not only for conventional ground-state calculations, but also for functionalities that go beyond conventional DFT calculations whereby explicit twoelectron Coulomb integrals and/or unoccupied KS states are needed [39][40][41][42]. In this regard, it is highly desirable to develop computational schemes within the NAO basis set framework that enables a convenient merging of firstprinciples and model-Hamiltonian based approaches. In case of DFT+U , several NAO-based implementations have been reported previously [43][44][45]. However, in case of DFT+DMFT, only one recent work within the pseudoatomic orbital basis sets (i.e., NAOs tailored for pseudopotential-based calculations) was reported by Sim and Han [46]. In their work, these authors proposed to use the so-called natural atomic orbitalsthe eigen-orbitals for the local density matrix -to define the local correlated orbitals. In the present work, we develop a projection scheme that allows to conveniently interface the NAO-based DFT codes and the DMFT impurity solvers and thus enables NAO-based DFT+DMFT calculations. We test the efficacy of such a scheme for the all-electron FHI-aims code [37] and the norm-conserving pseudopotential based ABACUS code [38], interfaced with three different DMFT impurity solver packages. Consistent results are obtained for transition metal compounds, including both correlated metals (SrVO 3 ) and Mott insulators (NiO and MnO), rare-earth systems (Ce metal and Ce 2 O 3 ), and actinides (Pu 2 O 3 and PuO 2 ). Our work thus extends the reach of the NAO-based numerical framework to treat strongly correlated materials. This paper is organized as follows. Sec. II focuses on the DFT+DMFT formalism. After introducing the general DFT+DMFT formulation, we present our definition of local correlated orbitals within the NAO framework. This is followed by a discussion of the self-consistency scheme of our DFT+DMFT implementation. In Sec. III, we benchmark the validity of our DFT+DMFT formalism and implementation over a wide range of strongly correlated materials, i.e., 3d -systems SrVO 3 , MnO and NiO, 4f -systems Ce metal and Ce 2 O 3 , and 5fsystems Pu 2 O 3 and PuO 2 . Sec. IV concludes this paper.

A. General DFT+DMFT formalism
Combining DFT with DMFT is not straightforward conceptually as the two theories are formulated in different forms and in different Hilbert spaces. In KS-DFT, one needs to self-consistently solve an effective singleparticle problem, whereby (for periodic systems) a KS Hamiltonian is constructed and solved separately at each individual k point in the Brillouin zone (BZ). On the other hand, DMFT is a real-space approach developed to solve the lattice Hamiltonian models in which the onsite two-electron interactions are explicitly included. Furthermore, as a first-principles approach, DFT accounts for all the chemical details in the systems and include all (at least valence) electrons in the calculations. In contrast, only the strongly correlated electrons sitting energetically in the vicinity of the Fermi level, relevant for the low-energy physics, are included in the model Hamiltonian and treated explicitly in DMFT. As such, key quantities like charge density, density matrix, and Green functions all have different representations in DFT and DMFT, impeding a straightforward combination of the two approaches. One solution to address this challenge, as detailed below, is to define a suitable projector to convert the key quantities back and forth between the two representations. As a side note, in the following we only focus on the general formalism of DFT+DMFT, without going into a detailed account of DMFT itself in this paper. Interested readers are referred to Refs. [7-10, 12, 47] for a detailed discussion of this theory.
Briefly, in the DFT+DMFT formalism, the Green functions defined in the correlated subspace are the central objects gluing the two theories. There are two kinds of interacting Green functions defined in the corrected subspace: One is the impurity Green function G imp determined by solving the impurity problem in Anderson impurity model, describing a single site with on-site Coulomb interactions and coupled to a meanfield electron bath. The another one, denoted as G loc , is the on-site term of the lattice Green function, which is essentially the the Green function of the correlated site described by the Hubbard model and is obtained by projecting the Green function in the full Hilbert space into the correlated subspace. The fundamental requirement in DMFT is that the "on-site" Green function of the lattice problem and that of the Anderson impurity model are equal, which is achieved via the DMFT self-consistency cycle.
In the DFT+DMFT scheme, the single-particle effective Hamiltonian in the full Hilbert space can be expressed asĤ whereĤ KS is the non-interacting KS Hamiltonian,Σ the self-energy which encodes all the many-body complexities arising from strong electron correlations and are nonzero only for correlated orbitals. Finally,Ĥ dc is the doublecounting term, which is introduced to discount the interactions among correlated electrons that have been included in a static mean-field level inĤ KS . Starting from the Hamiltonian Eq. (2), the key issue in the combination of DFT and DMFT is the downfolding (projecting) of the physical quantities from the full KS space to the local correlated space, and the upfolding (embedding) of these quantities from the local space to the full space. Mathematically, the local interacting Green function G loc can be obtained through a projection procedure. In the literature, two different projection procedures have been used. The first one can be seen as a "Hamiltonian projection", whereby a tight-binding Hamiltonian is first obtained from the KS Hamiltonian by projecting the latter into the correlated subspace, Here φ mk is the Bloch-summed local orbitals with m denoting an orbital index in the correlated subspace. In the second step, the local interacting Green function can be directly calculated from this tight-binding Hamiltonian as where N k is the is the number of k points in the first Brillouin zone, equivalent to the number of unit cells in the Born-von Kármán supercell. In Eq. (4), µ is the chemical potential, I an identity matrix, and [· · · ] −1 denotes the matrix inversion which is taken within the correlated subspace. Here and below, we employ the Matsubara Green function formalism where ω n = (2n + 1)π/β are the discrete frequency points along the imaginary frequency axis, with β = 1/k B T being the inverse temperature. The spectroscopy results on the real frequency axis are obtained via analytical continuations. The above procedure to compute the local interacting Green function is the so-called tight-binding Hamiltonian method.
In the second scheme, one first constructs the interacting Green function in the full Hilbert space as where the Hamiltonian matrix H ij (k) = χ i,k |Ĥ | χ j,k is the representation of the interacting Hamiltonian operator introduced in Eq. (2) within an arbitrary orthonormal basis set | χ j,k expanding the full Hilbert space. While the matrix elements of the KS Hamiltonian H KS within such a basis set are given straightforwardly, those ofΣ andĤ dc will be discussed later in Sec. II C.
Once the full-space lattice Green function is obtained, the local interacting Green function can be obtained through a projection procedure where the projector above is given by The construction of the local Green function via Eqs. (5)-(7) is known as the projector method. The two methods discussed above to construct the local Green function differ in that the tight-binding Hamiltonian method projects the Hamiltonian matrix while the projector method projects the Green function. In the projector method, the matrix inversion is carried out in the full Hilbert space, as indicated in Eq. (5), and thus the interactions between the correlated electrons and the rest electrons are retained to some extent. On the other hand, within the tight-binding method, if one wants to describe the interaction between correlated electrons and the rest ones, e.g., studying the charge transfer process between the correlated orbitals and the bath, one has to enlarge the Hilbert space of the tight-binding Hamiltonian to encompass the extra itinerant electrons, and it will inevitably increases the computational complexity. In this work, we choose the projector method in our DFT+DMFT implementation.
The other key quantity in DFT+DMFT calculations is the interacting impurity Green function G imp in Eq. (1), corresponding to the local propagator of the effective single-impurity Anderson model, which describes a single site coupled to a bath that mimics the lattice environment at a mean-field level. Formally G imp satisfies the following relationship where E imp is the impurity energy level, ∆ (iω n ) the socalled hybridization function characterizing the influence of the environment on the embedded impurity, and Σ imp (iω n ) the impurity self-energy. It should be noted that E imp and ∆ (iω n ) are matrices when there are multiple orbitals in the correlated subspace, which is typical in DFT+DMFT calculations. In this context, it is also customary to define a so-called Weiss Green function which acts as the dynamical (energy-dependent) mean field that the impurity electrons experience, and encodes essentially the same information as the hybridization function ∆ (iω n ). When the self-consistency in the DMFT loop is reached, the local Green function G loc mm (iω n ) and the local self-energy Σ mm (iω n ) in Eq. (4) and (6) are equal to the impurity Green function G imp mm (iω n ) and the impurity self-energy Σ imp mm (iω n ), respectively.
The Weiss Green function together with the local Coulomb interactions defines the Anderson impurity model, which can be expressed in terms of an action [7] where c † l (τ ), c n (τ ), etc., should be understood as the Grassmann variables, and U lmno is the on-site Coulomb interaction expressed within a set of local orbitals (labelled by l, m, n, o) spanning the correlated subspace. The action S is essentially the integration of the Lagrangian over the imaginary time. For given S, the impurity Green function can be calculated as where i runs over all m indices, and Z is the partition function The interacting impurity Green's function defined via Eqs. (10)(11)(12) can then be obtained through a variety of numerical approaches -usually termed as impurity solvers. Up to date, several types of impurity solvers have been developed, including the quantum Monte Carlo (QMC) [13,[48][49][50], non-crossing approximation (NCA) [51][52][53], one-crossing approximation (OCA) [54,55], exact diagonalization (ED) [56,57], numerical renormalization group (NRG) [58,59], etc. Among all these impurity solvers, continuous-time quantum Monte Carlo (CTQMC) [60][61][62][63] provides access to both high and low energy scales and is effective for a wide class of realistic material calculations. Nowadays, the CTQMC, especially the hybridization expansion based CTQMC (CT-HYB), is the most popular impurity solver employed in DFT+DMFT calculations. Once the interacting impurity Green function is determined, the impurity self-energy Σ imp (iω n ) can be ob- which is usually done within the impurity solvers. This impurity self-energy will be taken as the updated local self-energy and fed into Eq. (3) and (4) or (2) and (5), from which a new local Green function and consequently a new Weiss Green's function can be obtained. This is where a new iteration starts. This self-consistency loop keeps going until the self-energy reaches convergence or the local and impurity interacting Green function satisfies the self-consistency condition, i.e., Eq. (1).

B. Construction of the projector
Within the NAO basis set framework, it's natural to take the d or f -type NAOs that contribute most to the electronic states around the Fermi level as the the local correlated orbitals. These NAOs by construction are localized and atom-centered, and thus satisfy the usual requirement of correlated orbitals. In early DFT+DMFT implementations, analogous atomic-like orbitals like LM-TOs were used. LMTOs are minimal basis sets in the sense that for each angular momentum channel there is only one radial basis function. In contrast with LMTOs, the NAO basis sets are of multi-ζ character meaning that there are more than one radial functions per angular momentum, thus offering a more accurate description of the electronic structure. In the past, the DFT+U method has been successfully implemented within NAObased DFT codes [43][44][45], whereby it turns out to be a good practice to choose the most localized d or f basis functions as the correlated orbitals to apply the Hubbard U correction. Thus the most localized d or f orbitals seems to span the correlated subspace rather well. In practice, since NAOs centering on neighboring atoms are non-orthogonal to each other, certain orthorgonalization procedure is needed to generate an orthonormal local basis set, which is convenient for DMFT calculations.
Below we shall discuss our procedure to construct the projector and local correlated orbitals to facilitate DFT+DMFT calculations within the NAO basis set framework. In analogy to the transformation between Bloch orbitals and Wannier orbitals, we define the following Bloch-summed atomic orbitals as where φ I,m (r − τ I − R) is a NAO located at atomic site I in cell R. Here the magnetic quantum number m labels the different orbitals in the correlated angular moment channel, such as the five d orbitals or seven f orbitals. Since NAOs on neighboring atomic sites are nonorthogonal to each other, it is obvious that the Φ k I,m (r)'s defined in Eq. (13) with different m are also nonorthogonal. The key next step is to employ the Löwdin orthonormalization procedure to Φ k I,m , i.e., where is the overlap matrix. The newly obtainedΦ k I,m orbitals are then orthonormal by construction. Afterwards, a Fourier transform is applied to get the correlated orbital in real space, i.e., The orthonormality of W R I,m is also guaranteed by construction.
We choose the KS states |Ψ ik as the basis sets |χ ik (cf. Eq. (7)) to span the full Hilbert space, and then Eq. (6) becomes where with ik being the KS eigenvalues. The projector, Eq. (7), then becomes where the superscript 0 denotes the central unit cell. Since in DMFT calculations, only the "on-site" Green function, where m and m orbitals are located in the same unit cell is needed, the projector is designed to project the full Green function into the central unit cell, without losing generality. Formally, the projector for a given correlated atom I and a wave vector k is a fourth-order tensor, but since it is separable and symmetric, only a second-order tensor, i.e., the overlap matrix Ψ ik | W 0 I,m needs to be computed and stored in practical implementations.
The whole DFT+DMFT scheme requires the orthonormality of local orbitals representing the correlated subspace. In the language of the projector, it requires the projector to satisfy the following orthonormal condition In principle, this condition is automatically satisfied if the summation over i goes over all the KS bands. In practical DFT+DMFT implementations, however, one truncates the full KS states into a small subset around the Fermi level, which means i just runs over bands that are located in a chosen energy window around the Fermi level (in the following, these subsets of bands are denoted as C). This truncation destroys the completeness of |Ψ ik and thus the orthonormality of the projector. To deal with this issue, one can introduce an extra transformatioñ to orthonormalize the projector. The transformation matrix in the above equation is given bỹ which is nothing but the overlap between the projections of the orthonormalized local orbitals |W 0 I,m 's within the subspace C. Mathematically, the local correlated orbitals we use above to construct the projector can be explicitly expressed as In this form, our scheme is similar in spirit to the projective Wannier-orbital scheme proposed by Anisimov et al. [24] in the context of LDA+DMFT calculations. In our case, the most localized d or f NAO plays the role of the LMTO in the work of Ref. [24] There are a few advantages of using NAOs to define the local correlated space. Firstly, this choice is physically intuitive and technologically straightforward within the NAO basis set framework. We do not need to spend extra efforts to generate a set of local orbitals and make sure they are physically reasonable atomiclike orbitals. Secondly, from both the theoretical and technical perspectives, the choice of NAOs to define the local correlated space and the resulting projection scheme are suitable for all NAO-based packages. Especially, the key quantities that are required in this formalism, e.g., the KS wave functions and the overlap matrix of the basis functions, exist naturally in NAO-based DFT code packages, and hence no additional efforts are required to calculate these quantities. Thirdly, due to its high flexibility, our DFT+DMFT infrastructure can be interfaced with a new NAO-based DFT code without much effort. We hope it can serve as a platform to enable NAO-based DFT codes to do DFT+DMFT researches on strongly correlated materials. In this work, we implement the DFT+DMFT interface and test it with two NAO-based DFT codes using different techniques, i.e., the pseudopotential-based ABACUS code [38] and full-potential all-electron FHI-aims code [37].

C. DFT+DMFT self-consistency scheme
In this section, we will explain step by step our DFT+DMFT calculation procedure, according to the workflow outlined in the flow diagram depicted in Fig. 1.
Step 1. The DFT+DMFT calculation starts from wellconverged DFT band structures. To get high-quality KS orbitals |Ψ ik , a dense k-point mesh is usually needed.
Step 2. With |Ψ ik , the projector defined in Eq. (20) can be straightforwardly constructed. To this end, the most localized d or f orbital of the correlated atoms in the NAO basis set is used to construct |W 0 I,m .
Step 3. In the DFT+DMFT iteration loop, the frequency-dependent self-energy Σ(iω n ) is determined by the impurity solver at each iteration step. To start with, the initial self-energy is set to be equal to the doublecounting term, i.e.,Σ ij (k, iω n ) = 0. Here the following choice of the double-counting term is used. Here n I is the total number of correlated electrons associated with the correlated atom I and is fixed during the DMFT cycles. In the spirit of reducing the necessity of introducing additional empirical parameters, n I is given by projecting the KS orbitals in the subset C to the local subspace as where f ik is the occupation number of KS orbital Ψ ik . This double-counting scheme is similar to the so-called fixed double-counting [22,64] scheme, which is considered to be able to improve the stability of DFT+DMFT self-consistency loop [64] by fixing the value of n I . The difference is that the nominal number of strongly correlated electrons is specified manually in the fixed double-counting scheme, whereas in our case this number is determined using Eq. (24).
Step 4. Using the projector constructed in Step 2, we embed (upfold) the self-energy back to the selected KS space subset C, which is expressed as Step 5. During the DFT+DMFT self-consistency iteration, the electronic chemical potential needs to be adjusted according to the newly obtained self-energy at each iteration, to keep the number of electrons hosted by KS bands in C conserved. Within the DMFT cycle, this condition means that where β is again the inverse temperature 1/k B T . The summation of imaginary frequency ω n should run from negative infinity to positive infinity. However, in realistic calculations, to save computational cost, the explicit summation over Matsubara frequency points is only carried within a frequency window [−ω N , ω N ], where the contribution from the frequency points outside the window is treated approximately. This is enabled by making use of the asymptotic behavior of the self-energy, i.e., lim n→∞Σii (k, iω n ) =Σ ii (k, ∞), whereΣ ii (k, ∞) is a real value. Then the Eq. (27) is approximated by When the chosen cutoff Matsubara frequency ω N is high enough so that the asymptotic behavior of the self-energy is correct, this approximation is of good accuracy.
Step 6. As all information is secured, the local interacting Green function is constructed from Eq. (17) For some CT-HYB impurity solver, the imaginary time hybridization function is needed.
Step 7. Solve the impurity problem with the determined µ, E imp I , ∆ I (iω n ) and the given on-site Coulomb interaction through the impurity solver to obtain the new self-energy and impurity Green function. In this paper, we use the Kanamori [65,66] form Coulomb interaction in which only the density-density term is included, so there are no sign problem in our CTQMC calculations.
Step 8. Check whether the self-consistency is reached. If the self-consistency is reached, stop the DFT+DMFT calcution, else return to Step 4.

A. d -electron systems
We first benchmark our DFT+DMFT implementation on three prototypical strongly correlated d -electron systems -SrVO 3 , NiO, and MnO. For the DFT part, we carried out GGA calculations using two NAO-based code packages -FHI-aims [37] and ABACUS [38]. In FHI-aims, the default tight tier 1 basis set is used for V, Ni, Mn, O, and Sr atoms, and the corresponding cutoff radii of the basis functions are 6.0 A, 6.0 A, 6.0 A, 6.0 A, and 8.0 A, respectively. In ABACUS we use the SG15 optimized norm-conserving Vanderbilt (ONCV) multi-projector pseudo-potentials [67][68][69] and the corresponding optimized double-ζ plus polarization (DZP) atomic basis sets, which comprise 4s2p2d1f basis functions with a cutoff radius of 9.0 Bohr for transition metal atoms and 2s2p1d basis functions with cutoff radii of 7.0 Bohr for O atoms and 10.0 Bohr for Sr atoms. In all DFT calculations, the Perdew-Burke-Ernzerhof (PBE) GGA exchange-correlation functional was used [70]. As for the DMFT part, we employed three different CT-HYB impurity solvers to solve the single-site impurity problem. These are PACS [71] developed by one of the present authors Gang Li, iQIST developed by Huang et al. [72][73][74], and the one used in eDMFT [30] developed by Haule at Rutgers University [62,75], which is referred to as "Rutgers" below in this paper.

SrVO3
SrVO 3 has a simple cubic perovskite structure without magnetism. DFT fails to reproduce the upper and lower Hubbard bands observed in experiments, and this calls for advanced theoretical and computational techniques [8]. In the past, this material has been extensively studied both theoretically and experimentally [25,[76][77][78][79][80][81][82], which provides abundant reference results. Therefore, SrVO 3 is an ideal test example for DFT+DMFT implementations [24,26,27,46,66]. The V 4+ cation with only one 3d electron is located in the center of the octahedron formed by its six surrouding axial O 2− ligand anions. In the presence of an octahedral crystal field, the five degenerate d orbitals split into two subsets: Three-fold degenerate t 2g orbitals (i.e., d xy , d yz and d xz ), and two-fold degenerate e g orbitals (i.e., d z 2 and d x 2 −y 2 ). The only one 3d electron of V 4+ occupies the lower-energy t 2g orbitals with the higherenergy e g orbitals being empty. As a common practice, we only consider three degenerate t 2g orbitals in our DFT+DMFT calculation.
In DFT calculations, we use 11×11×11 k-point mesh generated by the Monkhorst-Pack method [83]. The Hubbard U and Hund J parameters are set to be 4.0 eV and 0.65 eV, respectively, following the choice in the literature [26,27,46]. The DMFT calculation is carried out at a temperature of 300 K. DFT calculations give a group of bands around the Fermi level with substantial 3d characterizes, which are well separated from other bands [26,27]. We enclose the six (twelve if spin degree of freedom is taken into account) bands crossing the Fermi level in the subset of KS bands C. One of the powerful strengths of Green-function based DFT+DMFT approach is that it can deliver physically meaningful single particle excitation spectral functions, as given by the the imaginary part of Green function. Such spectral functions can be directly measured by photoemission and/or inverse photoemission experiments. In Fig. 2, the spectral functions of 3d electrons of SrVO 3 calculated by combining two NAO-based DFT packages -FHI-aims and ABACUS, and three impurity solvers -PACS, iQIST and Rutgers are presented. Despite the considerable differences underlying the implementations of the DFT codes, as well as the DMFT impurity solvers, the calculated spectral functions are remarkably similar. Most importantly, they all successfully reproduce the typical three-peak structure: A quasi-particle band located at the Fermi level with the lower and upper Hubbard bands on the two sides, arising from the strong electron correlations which traditional LDA and GGA fail to capture. The main features of the spectral functions given by our DFT+DMFT calculations are in good agreement with previous theoretical results except for some small details such as the exact peak positions and intensities [24-27, 46, 66, 76, 78, 79]. Our results also correctly describe the strong particle-hole asymmetry of the lower and upper Hubbard band as the intensity of the occupied lower Hubbard band is much smaller than the empty upper Hubbard band. This is the consequence that the V 4+ , with only one 3d electron, is far away from half-filling.
Similar to previous DFT+DMFT studies [24-27, 46, 66, 76, 79], our spectral function results are in excellent agreement with the experimental photoemission spec-trum [78]. The first main feature lies in the reproduction of the lower Hubbard band at around −2.0 eV, which is the manifestation of the strong correlation among 3d electrons. The second feature, in quantitative agreement with experiment, is that the lower Hubbard band vanishes at ∼ −1.0 eV, and then the quasiparticle band starts to rise sharply.

MnO and NiO
Next we test our DFT+DMFT scheme on the late transition metal oxides MnO and NiO. Both MnO and NiO crystallize in face-centered cubic (FCC) NaCl-type structures where 3d orbitals are split into threefold degenerate t 2g and twofold degenerate e g orbitals as in the SrVO 3 case. In the MnO system, we adopt the parameters of U =5.0 eV, J =1.0 eV, and the temperature T  Hubbard model and ab-initio DFT, static DFT+U [6,20,21] and the more advanced dynamical DFT+DMFT [84,85] successfully open the gaps in those prototypical transitional metal oxides. Then obtained band gaps are in quantitative agreement with experiments, if proper parameters are used.
We then apply our DFT+DMFT implementation to both MnO and NiO, and the obtained spectral functions are presented in Fig. 3 and Fig. 4, respectively. Again in each figure, we present six panels containing results obtained by combining two NAO-based DFT codes and three DMFT impurity solvers. The most prominent feature in the DFT+DMFT spectra of both MnO and NiO is that a sizeable gap is opened up, arising apparently from the strong Coulomb interactions whose effect is properly captured by DMFT. One can further see from the fact that the number of 3d electrons of Mn given by n = − m G imp mm (β) is 5. Similar to the SrVO 3 case, the results presented in all six panels are very close, except that the peaks given by the Rutgers impurity solver is slightly more pronounced.
For NiO, the number of nominal 3d electrons is approximately 8, and DFT+DMFT calculations lead to a splitting of the 3d bands, with fully occupied t 2g orbitals and half-filled e g orbitals. Our results in general agree well with previous DFT+DMFT results reported in the literature [84,85]. However, Fig. 4 also reveals that both the DFT codes and the impurity solvers have certain influence on the obtained DFT+DMFT energy spectra. Although there are no qualitative differences, the shape and width of the left peak, and the depth of the dip between the first and middle peaks show noticeable quantitative differences. It is not entirely clear to us yet what factors caused such differences. We note that, however, NiO is a prototypical charge transfer insulator [86], and the competition between the strong Coulomb interactions among the Ni 3d electrons and the hybridization between the Ni 3d and O 2p electrons governs its underlying physics. Thus a complete DFT+DMFT treatment of NiO should also include O 2p electrons in the game. With the O 2p excluded, the DFT+DMFT calculations are probably more sensitive to the details in the band structures and numerical techniques behind the impurity solvers. Figure 3 and 4 only contain the spectral information of 3d electrons. Our DFT+DMFT implementation also allows for calculating the k-resolved spectral function in the KS orbital space, as determined by where the real-frequency self-energy is evaluated by analytical continuation of the imaginary frequency-selfenergy through the maximum entropy formalism in Ref. [87]. The corresponding results for MnO and NiO are presented in Figs. 5 and 6, respectively. Now for simplicity only results obtained using the two DFT codes combined with the Rutgers impurity solver are presented. The agreement between our theoretical gaps and experimental gaps [88][89][90] is rather satisfactory. The two k-resolved spectral functions of MnO is nearly iden-tical whereas again noticeable differences, in particular regarding the spectral weights of certain bands, exist for NiO. Nevertheless the main features, e.g., the energy positions and dispersion of the bands, are reasonably similar in Fig. 6(a) and 6(b).
NiO has been extensively studied experimentally, and reliable experimental photoemission spectra are available. In Fig. 7, we directly compare our DFT+DMFT spectrum of NiO with the experimental spectrum of Sawatzky and Allen [89] measured by x-rayphotoemission (XPS) and bremsstrahlung-isochromatspectroscopy (BIS) techniques. Again, the two sets of theoretical spectra are obtained using FHI-aims and ABACUS codes interfaced with the Rutgers impurity solver, respectively. For the spectrum below the Fermi level, corresponding to the (negative of) energy cost for removing a particle from the occupied levels, we reproduce the sharp peak around −2.0 eV and a local minimum at −3.0 eV given by XPS. For the spectrum above the Fermi level, corresponding to the (negative of) energy cost for adding a particle to the system, our curves match perfectly with the BIS result. Although small deviations of the shoulder peak at approximately −3.5 eV are visible, our DFT+DMFT spectrum, in general, agrees well with previous theoretical work [84] and the experiment data [89]. In ABACUS, the pseudopotential method is used and only the valence electrons are explicitly included. Furthermore, as given in the beginning of this section, the spatial cutoffs of NAOs used in ABACUS are much smaller than those used in FHI-aims. By contrast, in FHI-aims full potential is used with all electrons being included in DFT calculations. Although these two NAO schemes are rather different, we get nearly identical DFT+DMFT results of SrVO 3 and MnO and reasonably similar results for NiO. It proves the efficacy and robustness of our DFT+DMFT formalism and implementation. We expect that the DFT+DMFT formalism presented in this paper should work for other NAO-based DFT codes as well.

B. f -electron systems
The f -electron systems, including 4f lanthanides and 5f actinides, are an important class of strongly correlated materials characterized by partially filled f -type orbitals. These systems exhibit a variety of exotic phenomena, such as the heavy-fermion behavior, metal-insulator transition, Kondo physics, volume collapses accompanying phase transitions, etc. It has been generally accepted that traditional DFT calculations based on static meanfield-type approximations such as LDA and GGAs do not provide adequate accuracy for describing these physical scenarios. DFT+DMFT has proved to be a powerful approach to tackle these systems and achieved remarkable successes in the last two decades. In this section we test our DFT+DMFT implementation on several typical f -electron systems, including 4f systems like α and γ cerium (Ce) metal and Ce 2 O 3 , and 5f systems like PuO 2 and Pu 2 O 3 .
DFT calculations are carried out using FHI-aims with the default tight tier1 basis sets. For these basis sets, the cutoff radii of Ce, Pu, and O elements are 7 A, 6.0 A, and 6.0 A, respectively. Analogous to the case of 3delectron systems, the GGA-PBE functional [70] is used in the DFT calculations and 11×11×11 Monkhorst-Pack k-point mesh [83] is used for Brillouin-zone integration, which is deemed dense enough to obtain accurate DFT band structures. Unlike the practice in previous sections where results obtained using three DMFT impurity solvers are shown for comparison, here for simplicity only the results obtained using the Rutgers solver will be presented.

Ce metal
Despite its simple FCC structure and the fact that only one 4f valence electron is present per Ce atom, Ce exhibits spectacular physical properties, which has attracted considerable research interest over the last forty decades [66,[91][92][93][94][95][96][97][98][99]. At low temperature and ambient pressure, Ce crystalizes in the α phase (smaller volume), where the system is paramagnetic and shows Pauli-like magnetic susceptibility without forming local magnetic moments. At high temperature, Ce transforms into the γ phase (larger volume), which instead carries local magnetic moments and exhibits Curie-Weiss behavior of magnetic susceptibility. When increasing pressure or decreasing temperature, Ce undergoes the famous isostructural γ-α phase transition accompanied by a 14% volume collapse and a drastic change of magnetic properties [66,[93][94][95][96][97][98][99]. A theoretical understanding of this behavior poses a great challenge to condensed matter physics and has motivated lots of experimental and theoretical investigations. From the perspective of first-principles computations, standard DFT calculations are unable to give a proper description of this phase transition, and the approaches going beyond conventional DFT are required.
In this paper, the crystal structures of α and γ Ce phases are discerned by different volumes, i.e., 29 A 3 for α and 34 A 3 for γ phase, respectively, which covers the volume ranges during the α-γ transition [94]. The chosen energy window for KS subset C is [−5.0 eV, 5.0 eV]. The Coulomb interaction parameter U given by constrained DFT is 6.0 eV [66,98]. Here we used the value of U of 6.0 eV and Hund parameter J of 0.5 eV [96,98]. DFT+DMFT calculations are done at the temperature of 800 K.
Our DFT+DMFT spectra for Ce are shown in Fig. 8. First of all, the energy spectra of both phases show the characteristic three-peak structures of strongly correlated metals. One can also see that, going from the α phase to the γ phase, the central quasiparticle peak gets substantially suppressed, and its spectral weights transfer to the lower and upper Hubbard bands . Such results are in agreement with experimental observations [91,92] and previous DFT+DMFT calculations for Ce [94][95][96]98]. The relative strengths of the quasi-particle peak and the upper Hubbard peak of both α and γ phases are also captured qualitatively. However, a direct comparison of the calculated spectra with the experimental spectra reveals that the energy separation between the quasi-particle peak and upper Hubbard band is underestimated in our calculations, while that between the lower Hubbard band and quasi-particle peak is overestimated. Compared to previous DFT+DMFT calculations, the location of the lower Hubbard peak of our results is close to that reported in the literature [94,98] where the maximum of lower Hubbard peaks is around -4.0∼−3.0 eV. Concerning the energy positions of the upper Hubbard band, we found that the results reported in the literature also vary quite a bit, with those reported in Ref. [98] agreeing best with experiment. Our results are close to those reported in Ref. [96], where the upper Hubbard peaks are located at around 2.5 eV. In general, our DFT+DMFT spectrums are in quantitative agreement with experimental and previous theoretical results. The discrepancies between our results and experimental and previous theoretical results indicate that details of the DFT calculations, the definition of the projector, as well as the numerical implementation of the impurity solvers will still have appreciable influence on the outcomes of DFT+DMFT calculations. Further investigations along these lines are still needed. The resonant inverse photoemission spectroscopies (RIPES) are taken from Ref. [91], and the photoemission (PE) spectra are taken from Ref. [92].

Ce2O3
Ce 2 O 3 is a Mott insulator with a gap of about 2.4 eV [22,100] arising from the strong correlation among 4f electrons. This system has been calculated using a variety of approaches such as DFT+U, GW [101], and DFT+DMFT [22,28]. Ce 2 O 3 crystalizes in the hexagonal lattice with space group P3m1. The experimental lattice parameters a = b = 3.891 A, c = 6.059 A [102] are used in our calculations. The energy window for KS subset C is [−3.0 eV, 2.0 eV]. Constrained DFT calculation predicted the parameter U varying from 5.5 to 8.0 eV, and we adopt U as 6.5 eV and J as 0.5 eV. The temperature is set to be 300 K.
Here we calculate the total density of states (TDOS) of Ce 2 O 3 through where A(k, ω) is defined in Eq. (31). The result is shown in Fig. 9, our DFT+DMFT calculation gives a band gap of about 2.5 eV, which agrees well with the experimental and previous theoretical results [22,28,100,101]. For example, all predict that the lower Hubbard peaks are sharp and narrow whereas the upper Hubbard peaks are strong and wide. The interesting feature about the upper Hubbard part is that it emerges as a low plateau between 1 and 3 eV above the Fermi level, continued by a pronounced peak between 3 and 4 eV. Note that the peaks mainly composed of O p bands sitting below the lower Hubbard peaks does not show up here due to the fact that those O p characteristic bands are not included in the subset C in the DFT+DMFT calculations whereas they are retained in Refs. [22,28,101].

Pu2O3 and PuO2
Plutonium (Pu)-based oxides such as Pu 2 O 3 and PuO 2 are essential for the fuel components in current nuclear reactors, as well as transmutation of the minor actinides from spent nuclear fuels [103,104]. A clear understanding of the physio-chemical properties of plutoniumbased oxides is of key importance for the safe operation and development of nuclear reactor systems and nuclear waste reprocessing [105,106], and the correct description of the oxidation and reduction processes [107,108]. The physical, chemical, and thermodynamical properties of Pu-based oxides such as the chemical bonding and electronic structure are intimately related to the states of strongly correlated 5f electrons. Conversely, the relative tendency of delocalization versus localization of strongly correlated 5f electrons is extremely sensitive to the physical and chemical environment of the Pu atom [109]. The description of complicated behaviors of 5f electrons, e.g., whether they are settled in the delocalized or localized states, or in the intermediate regime, is out of the reach of standard DFT. In contrast, DFT+DMFT is becoming a promising approach that facilitates an indepth understanding of this type of materials.
In this paper, we look at two plutonium oxides, i.e., Pu 2 O 3 in its β phase [110] (simply called Pu 2 O 3 below) and PuO 2 [111] with our DFT+DMFT implementation. We use the Hubbard parameter U We calculated DFT+DMFT TDOS of Pu 2 O 3 through the method introduced in the case of Ce 2 O 3 . The TDOS here contains contributions from both the correlated orbitals and the rest (here spd) orbitals. The projected density of states (PDOS) belonging to correlated (5f ) orbitals, is evaluated by whereas the PDOS for the spd orbitals is obtained by taking the difference between the TDOS and ρ f (ω). The results are depicted in Fig. 10. Figure 10 indicates that the band gap of Pu 2 O 3 as determined by DFT+DMFT is about 1.7 eV, which is in good agreement with the experimental bands gap of 1.8 eV [28,112]. Additionally, the occupation analysis yields an average occupation number n f =5.0 for Pu 5f electrons through n = − m G imp mm (β), which is consistent with the chemical environment of Pu 3+ in Pu 2 O 3 .
In Fig. 11, the DFT+DMFT TDOS of PuO 2 is plotted. The bands gap is predicted to be about 2.5 eV, which is in good agreement with previous theoretical work [113] and the experimental band gap of 2.8 eV [114]. In PuO 2 , the peak of occupied Pu 5f electrons is noticeably sharper than that of Pu 2 O 3 , which indicates that the Pu 5f electrons in PuO 2 are likely more localized than the case of Pu 2 O 3 . This localization picture is in agreement with the bigger gaps in PuO 2 . The occupation analysis gives an average occupation number of n f =4.0 for Pu 5f electrons, which agrees with the chemical state of Pu 4+ in PuO 2 .

IV. SUMMARY
In summary, we developed and implemented a formalism that allows us to carry out DFT+DMFT calculations within the NAO basis set framework. For transition  metal compounds and f -electron systems, the most localized d or f -type NAO is used to define the local correlated subspace. Such a choice is physically intuitive and implementationally convenient for NAO basis sets.
Following what is usually done in the literature, only a subset of the KS bands C out of the full KS space around the Fermi energy, which hosts the majority of electrons with strong correlations, are enclosed in the definition of the projector. Our projector scheme is mathematically equivalent to the projective-Wannierfunction approach adopted in Ref. [24]. We implemented such a DFT+DMFT formalism by interfacing two NAObased DFT codes, i.e., FHI-aims [37] and ABACUS [38] with three DMFT quantum impurity solvers, i.e., PACS, iQIST [72][73][74], and the Rutgers impurity solver [30,62,75]. In particular, interfacing with the allelectron FHI-aims code allows one to study all types of strongly correlated materials over the periodic table.
Our DFT+DMFT formalism and implementation are testified for three typical series of strongly correlated materials, namely the 3d transitional metal compounds SrVO 3 , MnO and NiO, 4f materials including the Ce metal and Ce 2 O 3 , as well as the 5f actinides Pu 2 O 3 and PuO 2 . For SrVO 3 and MnO, our calculated one-electron removal and addition spectra are in good agreement with previously reported DFT+DMFT results and experimental data. Furthermore, the calculated results are rather robust against the use of different DFT codes or different impurity solvers. For NiO, the obtained DFT+DMFT results show noticeable dependence on the chosen DFT codes and/or impurity solvers, although the key spectroscopic features are captured by all calculations. For felectron systems, the characteristic three-peak structures are obtained for the Ce metal, whereas for the correlated insulators, the obtained DFT+DMFT band gaps are in good agreement with experiments. However, there remain issues calling for further investigations, like the energy separation of the quasiparticle peak and the upper Hubbard band for the Ce metal.
Admittedly, our scheme is still in its infancy, and at the quantitative level there are still issue to be sorted out. Further improvement would be necessary for more reliable descriptions of charge-transfer-type Mott insulators, and the intricate lanthanides and actinides. However, we consider that our attempt of developing an infrastructure that merge NAO-based DFT codes and DMFT-based techniques is rather rewarding. This will not only enable standard DFT+DMFT calculations within the NAO basis set framework, but also paves the way for developing more advanced schemes by combining beyond-DFT approaches like hybrid functionals [40,41,115] or GW [42], recently available in NAO-based DFT codes, with the DMFT machinery.