Ab Initio Full Cell GW+DMFT for Correlated Materials

Quantitative prediction of electronic properties in correlated materials requires simulations without empirical truncations and parameters. We present a method to achieve this goal through a new ab initio formulation of dynamical mean-field theory (DMFT). Instead of using small impurities defined in a low-energy subspace, which require complicated downfolded interactions which are often approximated, we describe a full cell $GW$+DMFT approach, where the impurities comprise all atoms in a unit cell or supercell of the crystal. Our formulation results in large impurity problems, which we treat here using an efficient coupled-cluster impurity solver that works on the real-frequency axis, combined with a one-shot $G_0W_0$ treatment of long-range interactions. We apply our full cell approach to bulk Si and two antiferromagnetic correlated insulators, NiO and $\alpha$-Fe$_2$O$_3$, with impurities containing up to 10 atoms and 124 orbitals. We find that spectral properties, magnetic moments, and two-particle spin correlation functions are obtained in good agreement with experiments. In addition, in the metal oxides, the balanced treatment of correlations involving all orbitals in the cell leads to new insights into the orbital character around the insulating gap.

DMFT, one must work within a diagrammatically clean formalism. In this context, it is natural to replace DFT with the GW approximation [18] as the low-level theory. The GW approximation (often employed in its one-shot form (G 0 W 0 )) [19,20] has been shown to fix many of the problems with semilocal density functionals (such as the underestimation of band gaps of weakly-correlated semiconductors) and thus appears a practical way to include the most important low-order long-range interaction diagrams. The combination with DMFT can then be formulated without double-counting by exactly subtracting the local GW contributions. The idea of self-consistently embedding the impurity self-energy and contributions to the polarization propagator arising from longrange interactions was proposed almost 17 years ago as the GW+(E)DMFT approximation [21,22], but only very recently have self-consistent implementations appeared [23,24]. However, while these developments are promising, applications have remained more limited than those with DFT+DMFT and have retained some problematic issues of that approach [25][26][27][28]. In particular, all current GW+(E)DMFT methods still require strongly downfolded interactions, because the impurity is restricted to the truncated low-energy subspace of a few correlated or orbitals ( Fig. 1(b)). Downfolding to a small number of strongly coupled orbitals is numerically challenging, and yields retarded interactions that either limit the applicable impurity solvers or which must be truncated or otherwise approximated. If one ignores the embedding of the polarization propagator to work purely with the bare interactions, one obtains the self-energy embedding theory (SEET) [29]. However, applications of this simpler approach in solids also remain limited [30]. Aside from these technical issues, in some more complex correlated materials, the local orbitals can be intertwined with other itinerant bands [31,32]. In such cases, even defining a set of local correlated orbitals can be difficult, and the quality of the calculation then depends sensitively on this choice [33].
A common origin of many of the above challenges is the definition of the impurity problem in terms of a small lowenergy subspace. This is done only to obtain as simple an impurity problem as possible, as motivated by model Hamiltonians, but it is not a requirement of the more general DMFT formalism. Consequently, in this work, we present a new formulation, which we term ab initio full cell GW+DMFT. In this approach we define the impurity to be the full unit cell -or even multiple unit cells of atoms -where each atom is described by a large localized set of atomic orbitals, covering the core, valence and high-energy virtual orbitals ( Fig. 1(a)). Since no low-energy subspace is identified, there is no downfolding, and we can simply use the full set of bare Coulomb interactions between the impurity orbitals, avoiding theoretical and numerical ambiguities. While conceptually simple, our full cell approach engenders two new complexities. The first is the need to set up the large impurity problem (for example to efficiently generate all the matrix elements) but this is enabled by technical advances we have made in the PySCF simulation platform [34] and our recently developed general ab initio quantum embedding framework [35,36]. The second is the need to solve the resulting impurity problem with a large number of orbitals. Here, the key insight is that many orbitals in the full cell impurity are only weakly correlated, and impurity solvers which take advantage of this, such as those used in molecular quantum chemistry [37,38] can then work very efficiently. In this work, we will use a coupled-cluster singles and doubles (CCSD) Green's function solver [39,40] carrying out self-consistency along the real-frequency axis. We apply the full cell GW+DMFT method to compute the spectral properties of Si as well as the spectra, magnetic moments and spin correlation functions of two correlated insulators, NiO and -Fe 2 O 3 , in their antiferromagnetic (AFM) phases. Our largest calculation in hematite uses an impurity of four Fe and six O atoms, giving rise to an unprecedentedly large ab initio DMFT impurity problem with 124 impurity orbitals.

II. THEORY AND IMPLEMENTATION
In the full cell GW+DMFT formulation, because the impurity cell contains all atoms in a crystal cell (or supercell), the effects of what would normally be thought of as long-range interactions on the polarization and self-energy from within the cell are all included. However, we will treat contributions from long-range interactions beyond the cell only at the level of the self-energy matrix of the crystal, computed at the one-shot G 0 W 0 level. Because of this, certain contributions to the polarization propagator involving interactions far from the cell, that would require the bosonic self-consistency of EDMFT [41], are omitted. Our formulation is designed to capture the main effects of polarization on the local strongly correlated degrees of freedom, while avoiding the full cost of the bosonic (polarization) self-consistency loop.
Given a periodic crystal, we start by performing a one-shot G 0 W 0 calculation on top of a mean-field reference (DFT or HF), using crystalline Gaussian atomic orbitals and Gaussian density fitting (GDF) integrals [42]. Because the G 0 W 0 approximation is reference dependent, we will denote the approximation G 0 W 0 @reference. The full G 0 W 0 self-energy matrix is computed in the mean-field molecular orbital (MO) basis along the imaginary-frequency axis [43,44]: where represents the 3-index electron repulsion integral (ERI) | | , is the Gaussian auxiliary basis, and and represent mean-field molecular orbitals (bands).
, and satisfy crystal momentum conservation: = − + , where is a lattice vector, and 0 ( , − ′ ) is the mean-field Green's function. The integration in Eq. 1 is carried out efficiently using a modified Gauss-Legendre grid [43] (100 grid points were used in this study). The polarization kernel ( , ′ ) is where and are occupied and virtual orbital energies respectively. Note that in a Gaussian basis formulation, the number of bands and size of auxiliary basis are significantly smaller than in plane-wave GW formulations [45], and because of this, the summations in Eqs. 1-2 run over all bands. To obtain the real-frequency G 0 W 0 self-energy, we perform analytic continuation. Here, we fit the self-energy matrix elements to -point Padé approximants ( = 18 in this work) using Thiele's reciprocal difference method [46]: To define the impurity problem, we first construct an orthogonal atom-centered local orbital (LO) basis. As in our previous work on ab initio HF+DMFT and density matrix embedding theory (DMET), we employ crystalline intrinsic atomic orbitals (IAOs) and projected atomic orbitals (PAOs) as the local orthogonal basis [35,36,47]. IAOs are a set of valence atomic-like orbitals that exactly span the occupied space of the mean-field calculations, whose construction only requires projecting the DFT/HF orbitals onto predefined valence (minimal) AOs. PAOs, on the other hand, provide the remaining high-energy virtual atomic-like orbitals that complete the atomic basis and capture the correlation and screening effects.
The impurity consist of all LOs (i.e. all IAOs and PAOs) within the impurity cell (crystal cell or supercell) with IAOs representing the core and valence orbitals and PAOs representing the high-energy virtual orbitals. This is illustrated in Fig. 1(a). The most expensive step in forming the impurity Hamiltonian is computing the bare Coulomb interaction matrix ( | ) for all orbitals within the impurity cell. However, using Gaussian density fitting, we can do this at relatively low cost (scaling asymptotically as ( 2 3 AO ), where , and AO are the numbers of points and auxiliary Gaussian and atomic orbitals within the impurity cell). We refer readers to Ref. [36] for a detailed algorithm. The impurity Hamiltonian (without bath orbitals) is thereforê with the one-particle interactioñ defined as the Hartree-Fock effective Hamiltonian with the double-counting term subtracted and imp is the impurity block of the mean-field density matrix. We then start the DMFT cycle with an initial guess of the impurity self-energy as the G 0 W 0 local self-energy: The G 0 W 0 local self-energy is computed in the LO basis within the impurity cell: and analytically continued to the real axis. Here, all local orbital indices ( , , , ) are within the impurity cell. The 3-index tensor is computed from a Cholesky decomposition of the impurity ERI: ( | ) = ∑ . Note that the polarization propagator is computed in the impurity orbital space, first in the imaginary time domain [48,49]: and then cosine transformed into imaginary frequency space. The hybridization self-energy is then computed: with the lattice Green's function defined as (9) and the full GW+DMFT self-energy defined as Here, is the chemical potential, which is adjusted during the DMFT self-consistency to ensure that the electron count of the impurity is correct. imp and ( ) are bare one-particle interactions for the impurity and whole solid.
In order to use a wavefunction (Hamiltonian-)based impurity solver, we discretize ( ). We discretize along the realfrequency axis [50] so that dynamical quantities (e.g., spectral functions) are obtained more accurately. To obtain the discretization, we optimize bath couplings { ( ) } and energies { } to minimize a cost function over a range of real-frequency points: where is the number of bath energies and is the number of bath orbitals per bath energy, and we use a broadening factor = 0.1 a.u. . The bath degrees of freedom are truncated by only coupling bath orbitals to the valence IAOs, further reducing computational and optimization costs. The full embedding problem with both impurity and bath orbitals is thus defined from the Hamiltonian We solve for the ground-state and Green's functions of the impurity Hamiltonian using a CCSD Green's function solver. At the singles and doubles level, CC may be viewed as generating ring, ladder, and coupled ring-ladder diagrams, and is known to be accurate in a variety of settings, including simple metallic and ordered magnetic states in ab initio calculations [51][52][53], and across weak to strong couplings when employed with small cluster DMFT impurities in Hubbard-like models [39]. (Note that the CCSD solver is certainly not the optimal solver for all problems, and other solvers will be explored in future work). From the CC Green's functions, we obtain an updated impurity self-energy imp ( ), and from this the DMFT cycle (Eqs. 8-12) is iterated until convergence between the impurity and lattice Green's functions:

III. RESULTS
We first apply our method to crystalline silicon. Although Si is considered a weakly-correlated semiconductor, it is still a challenging system for many DFT functionals (such as LDA and GGA) which do not yield accurate band gaps [54]. Oneshot G 0 W 0 on top of LDA or GGA is known to significantly improve the band structure, although this relies somewhat on the cancellation of errors [55]. Such a small band-gap system also poses challenges to quantum embedding methods that start from a local correlation picture, such as DMFT [56], due to the long-range nature of its statically screened Coulomb interaction, which must be included in the treatment. The full cell GW+DMFT results for Si are presented in Fig. 2. We used the GTH-PADE pseudopotential [57] and GTH-TZVP basis [58], and a 6 × 6 × 6 Γ-centered -point sampling. The impurity was defined as the unit cell of 2 Si atoms with 34 local orbitals (2 2 3 3 3 4 4 for Si), and 128 bath orbitals were used. As known from other G 0 W 0 calculations [59] and as seen in Figs. 2(a) and 2(b), the meanfield starting point strongly affects the quality of the G 0 W 0 results; G 0 W 0 @PBE gives an accurate band gap of 1.09 eV when compared to the experimental value of 1.17 eV [60], while G 0 W 0 @HF overestimates the band gap, giving 2.04 eV. GW+DMFT predicts the band gap of Si to be 1.01 eV (@PBE) and 1.39 eV (@HF), largely removing the reference dependence of G 0 W 0 , due to the more complete inclusion of diagrams from interactions within the unit cell. The spectral function is also greatly improved in GW+DMFT compared to G 0 W 0 @HF. In earlier HF+DMFT calculations [36] (see inset of Fig. 2(b)), we found the band gap to be too large by 0.5 eV, and this quantifies the effect of the long-range correlations in G 0 W 0 on the band gap of Si. From Fig. 2(c), we note that GW+DMFT@PBE maintains the accurate band structure of G 0 W 0 @PBE, in contrast to self-consistent GW, which is known to lead to worse results than G 0 W 0 itself [61,62]. We next show the results of full cell GW+DMFT in Fig. 3 for a strongly-correlated insulator, NiO, in the antiferromagnetic phase. The GTH-PADE pseudopotential and GTH-DZVP-MOLOPT-SR basis set [63] were used with a 6 × 6 × 6 Γ-centered -point sampling defined with respect to the antiferromagnetic cell (2 NiO units). (As an estimate of the remaining finite size error, the difference between the G 0 W 0 @PBE gaps for 4 × 4 × 4 and 6 × 6 × 6 -meshes is only 0.1 eV). We used the antiferromagnetic cell of 2 NiO units along the [111] direction as the impurity, corresponding to 78 impurity orbitals (3 3 3 4 4 4 4 for Ni and 2 2 3 3 3 for O) and 72 bath orbitals in the DMFT impurity problem. As seen in Fig. 3(a), G 0 W 0 @PBE severely underestimates the band gap at 1.9 eV, even when using a spin symmetry broken PBE reference. Meanwhile, the valence spectrum of G 0 W 0 @PBE does not agree well with the experimental photoemission spectrum [64]. GW+DMFT, on the other hand, predicts a band gap of 4.0 eV and a magnetic moment of 1. 69 , both in very good agreement with the experimental values of 4.3 eV [64] and 1.77-1.90 [65,66]. More interestingly, our GW+DMFT DOS captures the experimental twopeak structure of the valence spectrum around −2 and −3 eV. A detailed analysis of the spin-orbital-resolved local DOS in Fig. 3(b) reveals that this two-peak structure results from the splitting of the majority and minority spin components of the Ni-2 orbitals, and is a signature of the AFM phase, as it does not arise within the paramagnetic phase [67]. From the local DOS, we can also conclude that NiO is an insulator with mixed charge-transfer and Mott character, with a valence band with contributions from Ni-2 , Ni-and O-2 , and a conduction band that is mainly of Ni-character.
In Figs. 3(c) and 3(d), we further analyze the character of the conduction band minimum (CBM) and valence band maximum (VBM) in the Brillouin zone using the momentumresolved DOS. We find that the lowest conduction band has strong Ni-4 and O-2 character at the Γ point (CBM), which was not discussed in many earlier DMFT calculations [28,68,69] which focused on the Ni-3 and O-2 orbitals and thus did not include Ni-4 (or O-2 ) orbitals in the impurity (although see Ref. [70] for a notable exception), unlike our full cell GW+DMFT treatment. At the Z point (VBM), we find that the highest valence band has significant O-2 and Nicontributions, with very little Ni-2 character. This is very different from the local DOS, where the Ni-2 has dominant weight in the valence bands. We confirm this by plotting the spatially-resolved DOS of NiO in the (001) plane in Fig. 4. We see that at the first valence peak and around the Ni atoms, the local spatial DOS has a Ni-2 ( ) orbital shape, while the momentum-resolved spatial DOS (at the Z point) has a Ni-( 2 − 2 ) orbital shape. Further, the weight of the DOS around the O atoms in Fig. 4(b) is considerably larger than in Fig. 4(a). Since our impurity includes two NiO units, we can also look at correlations across the cells. We computed the spin-spin correlation function for the Ni atoms within the impurity problem: ∑ We found ⟨ ⋅ ⟩ between two Ni atoms to be -0.707. Both ⟨ ⟩ and ⟨ ⟩ contribute almost zero spin correlation, and the uncorrelated value ⟨ ⟩⟨ ⟩ is -0.710, suggesting that quantum spin correlations are weak and NiO is close to a classical Ising magnet. This is consistent with experimental measurements of the critical behavior of the magnetic phase transition in NiO [71] and our previous ab initio DMET study [35].
We next turn to study a second strongly-correlated insulator, hematite ( -Fe 2 O 3 ), in the AFM phase. We take the impurity to be the complete AFM unit cell, including 2 Fe 2 O 3 units (Fig. 1), with a "+ − +−" type AFM ordering of the Fe spins. Because of the large impurity size, we used a DZVquality basis (GTH-DZV-MOLOPT-SR, 3 3 3 4 4 4 5 for Fe, 2 2 3 3 for O), leading to an impurity problem with 124 impurity and 48 bath orbitals. The small number of bath orbitals is due to the current numerical limitations of our CCSD solver. However, since our bath orbitals are only coupled to the valence impurity orbitals, and we aim to reproduce the hybridization only in a window near the Fermi surface (±0.4 a.u.), the bath discretization error is not too severe. The 3 3 orbitals of Fe were treated as frozen core orbitals (i.e., uncorrelated) in the CCSD solver. The GTH-PBE pseudopotential and 4 × 4 × 4 Γ-centered -point sampling were employed. As presented in Fig. 5(a), G 0 W 0 @PBE severely underestimates the band gap at 0.5 eV, compared to the experimental value of 2.6 eV [72]. G 0 W 0 with the hybrid functional PBE0 slightly overestimates the gap (3.4 eV), but the spectrum does not agree well with experiment, and in particular, the features of the G 0 W 0 @PBE0 DOS are too sharp around -7 and 3.5 eV (Fig. 5(b)). GW+DMFT improves the G 0 W 0 @PBE spectrum significantly, especially in the valence region, although the band gap (1.5 eV) is still too small. From the orbital-resolved DOS in Fig. 5(c), we find that the main improvement comes from the spectral positions of the majority spin component of the Fe-3 orbitals and O-2 orbitals. G 0 W 0 @PBE mistakenly predicts the Fe-3 valence spectrum to lie close to the Fermi surface and that Fe 2 O 3 has considerable Mott insulating character. However, GW+DMFT shifts the majority-spin Fe-3 DOS to lower energies, consistent with previous DFT+DMFT calculations [73,74]. Because of this correction, GW+DMFT obtains a more accurate Fe magnetic moment than PBE (4.23 compared to 3.71 with the experimental moment being 4.64 [75]). We find that the valence band spectrum is dominated by O-2 near the Fermi surface, indicating that Fe 2 O 3 is in fact a pure charge-transfer insulator, with almost no Mott insulating character. This is in contrast to DFT+DMFT calculations [73,74] that find a sizable Fe-3 contribution to the valence band maximum. We attribute this disagreement to the full cell GW+DMFT treatment where both O-2 orbitals and Fe-3 are treated on an equal footing at the impurity level, which thus allows for a more accurate balancing of their relative contributions to the spectral weight.
Starting from a PBE0 reference, GW+DMFT finds a slightly larger band gap (3.9 eV) and magnetic moment (4.37 ) than G 0 W 0 (3.4 eV and 4.20 ). The overly sharp peaks of the G 0 W 0 @PBE0 spectrum around −7 and 3.5 eV are corrected by GW+DMFT, which broadens the Fe-3 peaks as shown in Fig. 5(d). In summary, it appears we achieve a good description of the photoemission spectrum for Fe 2 O 3 within the full cell GW+DMFT, although a fully quantitative prediction of the band gap is not attained. Given that G 0 W 0 only provides a minor correction to the underlying DFT band gap in this system, the likely culprit is the insufficiency of the G 0 W 0 approximation in describing the long-range interactions in Fe 2 O 3 .

IV. CONCLUSIONS
In this work, we introduced a full cell GW+DMFT formulation for the ab initio simulation of correlated materials. The primary strength of this approach is that it entirely avoids the problem of selecting a low-energy subspace, and consequently the uncontrolled errors introduced either by downfolding the effective interactions within the subspace, or via DFT double counting. The resulting method is then both parameter free and can easily treat all interactions. We showed that full cell GW+DMFT can be applied to systems using impurity cells of up to 10 atoms in calculations of the spectral properties of Si, NiO and -Fe 2 O 3 , obtaining for most quantities, results of good quantitative accuracy. By defining the impurity to comprise all orbitals in the AFM supercells of NiO and -Fe 2 O 3 , we also showed how the full cell approach can cleanly differentiate between different amounts of charge-transfer and Mott insulating character and the orbital character around the gap, as both metal and non-metal orbitals enter into the impurity problem on an equal footing. Overall, our calculations demonstrate the potential of the full cell GW+DMFT approach for studies of more complicated materials, en route towards a fully predictive theory of correlated materials.