Parameter-free hybrid functional based on an extended Hubbard model: DFT+U+V

In this article, we propose an energy functional at the level of DFT+U+V that allows us to compute self-consistently the values of the on-site interaction, Hubbard U and Hund J, as well as the intersite interaction V. This functional extends the previously proposed ACBN0 functional [Phys. Rev. X 5, 011006 (2015)]. We show that this ab initio and self-consistent pseudo-hybrid functional yield improved electronic properties for a wide range of materials, ranging from sp materials to strongly-correlated materials. This functional can also be seen as an alternative general and systematic way to construct parameter-free hybrid functionals, based on the extended Hubbard model and a selected set of Coulomb integrals, and might be use to propose novel approximations. By extending the DFT+U method to materials where strong local and nonlocal interactions are relevant, this work opens the door to the ab initio study the electronic ionic and optical properties of a larger class of strongly correlated materials in and out of equilibrium.


I. INTRODUCTION
During the last few decades, density functional theory (DFT) has emerged as one of the reliable and efficient numerical method to simulate a wide range of materials. However, it is well known that the most employed local and semilocal functionals suffer from many problems, in particular the so-called "delocalization problem", which prevent using DFT with usual functionals on materials where strong local electron-electron are taking place 1,2 . More advanced functionals, such as some based on the meta-generalized gradient approximation, or hybrid functionals can solve some of these problems, but are still not ideal for strongly-correlated systems. In order to overcome this problem, one well established method is to perform dynamical mean field theory (DMFT) calculations 3,4 , possibly using the outcome of a DFT calculation, a method usually referred as DFT+DMFT, which has become the state of the art method to treat strongly correlated materials 5 . In this framework, the effective electronic parameter describing the local interaction, the Hubbard U , is computed within the framework of constrained random-phase approximation (cRPA) [6][7][8][9] .
As an alternative effective approach, DFT+U method originally proposed by V. Anisimov, A. Lichestein, and coworkers 1,10-12 has become a well established and successful way to improve the treatment of correlated solids upon DFT, without the numerical burden of DFT+DMFT or GW+DMFT. In order to correct the over-delocalization of the electrons, it was proposed to include an energy penalty U for the localized 3d or 4f orbitals orbitals, in the spirit of the mean-field Hubbard model. 1,2,[10][11][12] The success of the DFT+U method mainly originates from the simplicity of the method, its relative low computational cost, and the fact that it can predict the proper magnetic ground state of chargetransfer and Mott insulators. 1 Of course, due to the strong simplifications involved in the DFT+U method, it has some intrinsic deficiencies, such as yielding infinite life-times for quasi-particles, or opening the bandgap only by making a long-range order in magnetic materials. The DFT+U approach is also not applicable to some really strongly-correlated systems and for these systems once needs to go beyond DFT+U and use DMFT or cluster DMFT frameworks. In this work, we focus on materials with not so strong correlations, but correlated enough such that the standard functionals do not capture the right correlation effects. The DFT+U approach still remains very attractive when it comes to the calculation of larger systems, such as twisted bilayer systems with small twist angles 13 , or for out-of-equilibrium situations, using realtime TDDFT+U 14,15 . It also improves the description of the optical properties of some correlated materials within linear response 16 .
Here, we develop an efficient numerical approach tailored towards materials for which not only local correlations are important, but also nonlocal correlations, which means a strong interaction between neighboring localized electrons. One example of such systems are charge-ordering insulators, such as Fe 3 O 4 17 , for which an electron is delocalized on two sites and hence the Mott-Hubbard localization cannot occur. Nonlocal interaction also play a key role in low-dimensional systems, such as ad-atoms on Si(111) surfaces 18 or sp electron systems such as graphite and graphene [19][20][21] . Finally, the role of the intersite interaction is strongly debated for high-T c superconductors and dictate many properties 22 . In all these systems, the intersite interaction plays a decisive role and the related self-interaction error contained in (semi-)local functional of DFT is crucially hampering the capability of these functionals. The low-energy physics of these systems is well described by an extended Hubbard model. In particular, if we only account for charge interaction between neighboring sites, the extended Hubbard Hamiltonian reads where σ denotes the spin index, t ij are the hopping matrix elements, U represents the on-site interaction, and V ij are the nonlocal Coulomb matrix elements between the neighboring sites i and j.
Our goal is to derive the expression of an energy functional containing at the same time U and J describing the multi-band on-site interaction, and the intersite interaction V describing the charge interaction between the different atomic sites. The functional can be seen as an hybrid functional, in which the Kohn-Sham orbitals are expended into a basis of atomic orbitals, including the onsite terms, and some of the intersite terms. At variance with most of the proposed hybrid functionals, we do not use a mixing parameter to determine the weight of the exchange interaction, but we base our approach on the extended Hubbard model and the related DFT+U + V scheme, which allows us to propose a fully ab initio and self-consistent scheme, which yields an estimate of the U , J, and V effective electronic parameters. Moreover, we use an approximate double-counting term, as commonly done for DFT+U , which does not requires any parameter to be adjusted. The self-consistency in the calculation of the Hubbard U can be crucial in the case of transition-metal complexes 23 , and we expect this to be equally true for the intersite interaction. In our approach, the on-site Hubbard U , Hund J and the intersite interaction V are all evaluated at the same time, ab initio and self-consistently, to ensure the consistency of our approach. It is important to stress that except for the additional cost of computing more Coulomb integrals at the beginning of the calculation, our approach does not represent a major extra cost compare to the usual DFT+U and the ACBN0 functional 16,24 . This makes our method very attractive. Another important application of method is the possibility to directly extend it to the time-dependent case, only assuming the adiabatic approximation, or to couple to other degrees of freedom, such as phonons.
This paper is organized as follow. First we briefly review the DFT+U and the DFT+U + V methods in Sec. II. Then we present our generalization of the ACBN0 functional, the extended ACBN0 functional, in Sec. III. We then test our new functional on different system and compare our results with prior works. Finally, we draw our conclusions in Sec.V.

II. DFT+U AND DFT+U
where E ee is the electron-electron interaction energy, and E dc accounts for the double counting of the electronelectron interaction already present in E DFT . Although an exact form of the double-counting term was recently proposed in the context of DFT+DMFT, 25 this doublecounting term is not known in the general case and several approximated forms have been proposed along the years. The E ee and E dc energies depend on the density matrix of a localized orbitals basis set {φ I,n,l m }, which are localized orbitals attached to the site I. In the following we refer to the elements of the density matrix of the localized basis as occupation matrices, and we denote them {n I,σ mm }. Combining these two expressions, we obtain the E U energy to be added to the DFT total energy, which only depends on an effective Hubbard U parameter where I is an atom index, n, l and m refer to the principal, azimuthal, and angular quantum numbers, respectively, and σ is the spin index. In the case of a periodic system, the occupation matrices n I,n,l,σ mm are given by where w k is the k-point weight and f σ nk is the occupation of the Bloch state |ψ σ n,k . Here, |φ I,n,l m are the localized orbitals that form the basis used to describe electron localization andP I,n,l mm is the projector associated with these orbitals, usually definedP I,n,l mm = |φ I,n,l m φ I,n,l m |. The definition of the orbitals and of the projector will be discussed in more details below. In the following, we omit the principal quantum number n for conciseness.
Recently an extension of the DFT+U method was proposed by V. Leiria Campo Jr and M. Cococcioni 27 , in order to account for the intersite electronic interaction V , in the spirit of the extended Hubbard Hamiltonian 22 , accounting only for the charge interaction between neighboring sites. In a similar spirit, Belozerov and coworkers proposed a LDA+DMFT+V approach that they applied to the monoclinic phase of VO 2 28 . The Hubbard U is usually defined as the averaged of the on-site interactions of the localized orbitals where V ee is the screened Coulomb interaction. In a similar way, the most akin definition of the averaged inter-site interaction V is defined as 27 (6) For conciseness, we omit below the quantum numbers in our notation and refer in the following to V I,l I ,l as V IJ . The expression for E ee and E dc become for DFT+U +V 27 where * IJ denotes that for each atom I the sum runs over its neighboring atoms J. This definition uses a generalization of the occupation matrix where it is clear that n IIσ mm is the usual occupation matrix n Iσ mm as defined in Eq. 4. Combining the previous expressions, we arrive to the energy E U V that must be added to the DFT energy This energy is the expression for the DFT+U + V proposed in Ref. 27, and is invariant under rotation of the orbitals of the same atomic site. This is a generalization of the work of Dudarev et al. 26 , where the double-counting expression is a generalization of the fully-localized limit (FLL) double-counting of DFT+U . The motivation of this specific expression for the intersite interaction was done in Ref. 27 is therefore not discussed here. Below we show how it is possible to extend the work of Agapito et al. 24 to evaluate the average intersite interaction V ab initio and self-consistently, in a form of a pseudo-hybrid calculation. 24, an approximation to the electron interaction energy named ACBN0 functional was proposed, allowing for an efficient ab initio evaluation of the DFT+U energy, which can be seen as a screened Hartree-Fock evaluation of the on-site U , or equally as a psuedo-hybrid functional in which the (screened) Hartree-Fock energy is included only on a selected localized subspace. We propose here an extension of this approach, which includes not only the on-site interaction, but also the charge exchange between two sites. Below, we refer to this new functional as the extended ACBN0 functional.
In our generalized functional, the electron interaction energy is given where we restricted only to the charge interaction between two sites, neglecting other two-site interaction and also three-and four-site interactions, as this is most likely to be the largest contribution to the energy 27 . In Eq. (9), the renormalized occupation matricesn I,σ mm and occupationsN I,σ ψ nk are respectively given bȳ where the sums run over all orbitals of the system owning the quantum numbers n and l, and being attached to atoms of the same type as the atom I, as this quantity is similar to the Mulliken charge of atom I. 24 Here we defined a generalized renormalized occupation matrices where we introduced an ansatz similar to the one of the original paper of the ACBN0 functional namely a degree of screening in thisn IJ,σ mm . It is clear that for the on-site case, this expression reduces to the expression of Ref. 24. Let us comment on the motivation for this renormalization of the generalized occupation matrix. In the case of on-site interaction, the original motivation was that if a wavefunction is fully delocalized, i.e. not occupying the localized subspace, it should gives a U of zero, and the ACBN0 functional should reduce to the (semi-)local functional used to describe the itinerant electrons. By including the weighting coefficient, which is the Mulliken charge of the set of orbitals, the effect is enhanced 24 . Here, as we consider an intersite interaction, the idea is the opposite. In the atomic limit, for which the wavefunctions are not delocalized over different sites, there should be no intersite interaction. This is well seen by the fact that for the Coulomb integral of the form (φ I m φ I m |φ J m φ J m ) to be nonzero, there must be an overlap of the charge density originating from the two atomic sites. Hence, if a wavefunc-tion is not delocalized over the two considered sites, the productN Iσ ψ nkN Jσ ψ nk vanishes, and the wavefunction does not contribute to the intersite V . Importantly, thanks to this weighting coefficient, if a wavefunction is fully localized on one site, the wavefunction does not screen the interaction. In the context of cRPA, the screening is given by the "rest" of the system, i.e., the electrons in the localized subspace do not contribute to the screening. This property is preserved by the renormalization factor of the ACBN0 functional, as well as in our extended ACBN0 functional, which at least partly explains the success of this approach.
From Eq. (9), the effective inter-siteV IJ is given, for J a neighbor of the atom I, bȳ This expression is the main result of this section. With it, one can evaluate the intersite interaction ab initio and self-consistently, similarly to what is done for the effective U in the ACBN0 functional. We also derived the expression of the extended ACBN0 functional for the case of noncollinear spins, as presented in Appendix A.
So far, we remained elusive on the orbitals used to define the localized subspace. In the Octopus code 16 , we construct the localized orbitals {φ I,n,l m } by taking the radial part of the pseudo-atomic wavefunctions given by the pseudopotential files, and multiplying it by the usual spherical harmonics, in order to obtain the pseudoatomic orbitals. More precisely, in case of periodic solids, we use in all the above equations not the isolated localized orbital, but the Bloch sums of the localized orbitals, which read as which amounts for introducing a phase factor when the sphere on which the atomic orbital is defined crosses the border of the real-space simulation box. Here N corresponds to the number of unit cells forming the periodic crystal, i.e., the number of k-points of the simulation. The projection to Kohn-Sham states selects a single momentum, explaining why the k-point index was not specified in the above equations. Also note that the normalization factor in the Bloch sum vanishes when the we use the periodicity of the crystal in computing the sum over the entire crystal, which reduces to a sum over the unit cell without the normalization factor.
In order to be able to treat various type of solids, including weakly correlated solids such as Si, we also implemented the Löwdin orthonormalization procedure, which transform the set of non-orthogonal localized orbitals {φ I,n,l m } into an orthonormal set of localized orbitals where i and j indices run over all the considered orbitals, and the overlap matrix for the set of considered orbitals is (S k ) ij = φ ik |φ jk . Importantly, using these orthogonalized orbitals, we obtain a trace of the on-site occupation matrix is consistent with the Mulliken population analysis and leads to the exact same trace of the occupation matrix than the dual projector defined in Ref. 29 for non-orthogonal basis set. Due to the periodicity of the Bloch sums of the localized orbitals, we only need to compute projections on the orbitals of the atoms inside the simulation box, irrespective of the number of neighboring atoms considered. This is also the case for a simpler DFT+U calculation, making the cost of DFT+U + V calculation only mildly more expensive to the one of a more standard DFT+U calculation, at variance with the original formulation of DFT+U + V , which requires the construction of larger supercells to include further neighbors 27 . We stress here that the orthonormal orbitals given by Eq. 14 are not the one we use for computing the Coulomb integrals, as they are periodic orbitals, and Coulomb integrals computed using these orbitals would contain both on-site and intersite overlaps, which is not what we want here. For this reason, the Coulomb integrals are computed before performing the orthonormalization procedure, from the atomic orbitals of the pseudopotential. More precisely, each of them are evaluated on a portion of the grid, and the Coulomb integrals are computed on the union of these two spheres, using a non-periodic Poisson solver. This is sketched in Fig. 1, in which the violet points correspond to the grid points obtained from the union of the two spherical meshes centered on two atoms (indicated by green crosses). This is the points we used in our implementation. Another possible choice would be to use a single large sphere centered at the middle of the two atoms. This is indicated in red in Fig. 1. This choice obviously leads to more grid points and we found that there is no difference in between the two choices, as the atomic orbitals rapidly decay away from the center of the atom. Finally we note that a formulation in terms of Wannier functions would look very much the same as the one we have presented here.

A. ZrSiSe
The nodal-line semimetal have received a lot of attention recently. Due to the vanishing density of states at the Dirac or Weyl points, the screening of the Coulomb interaction is altered, and long-range Coulomb interaction is a crucial ingredient in the description of these materials 30 Hence, the nodal-line semimetal exhibit strong nonlocal correlations 31,32 , which are not captured by a local Hubbard U as used in DFT+U . This represent therefore an interesting potential application of our extended ACBN0 functional. In order to benchmark our functional, we decided to investigate ZrSiSe, and use hybrid functional as a reference to compare with. In Ref. 31 and 33, it was shown that hybrid functional calculations with a fraction of exact exchange of 7% reproduce best the experimental results, compared to the standard fraction of 25% used in the HSE06 functional 34 . We employed the experimental lattice constant of 3.623Å and we sampled the real-space using a spacing of 0.3 bohr and the Brillouin zone using a 7 × 7 × 3 k-point grid. We considered localization of the d orbitals of Zr and included the interaction with the first nearest neighbors, and we employed the localdensity approximation for the DFT exchange-correlation  Fig. 2. We found that DFT+U + V has almost the same band dispersion than the one obtained from the hybrid functional close to the Weyl point, showing the validity of our functional in describing this material. We recently applied our functional to reproduce the measured angle-resolved photoemission spectrum of ZrSiSe 33 including spin-orbit coupling. Formulae for noncollinear spins are presented in Appendix A

B. Graphene and graphite
Low-dimensional sp materials received a lot of attention recently as they display strong local and nonlocal Coulomb interaction 20,21 . In Ref. 20, authors reported cRPA calculations of the on-site and intersite interaction for graphene at half-filling. Similar calculations were performed in graphene and graphite in Ref. 21. In order to illustrate the flexibility and robustness of our implementation, we compute the intersite interactions for graphene, treated here with mixed periodic boundary conditions, and graphite, which is a fully periodic material. We use a 15 × 15 k-point grid to sample the two-dimensional Brillouin zone of graphene and a 12 × 12 × 4 grid in the case of graphite. We employ an in-plane lattice constant of 2.47Å and an out-of-plane constant of 6.708Å for graphite. The real-space grid is sampled by a spacing of 0.45 Bohr and we employed the norm-conserving Pseudodojo pseudopotential 35 . We compare the previously reported values to the ones obtained by our functional in Tab. I. A major difference is that the in the present calculations all p orbitals are considered, whereas prior studies only considered p z orbitals. As a result, both onsite and intersite interaction are found to be smaller in our case that in previous works 20,21 , which is fully compatible with considering more orbitals in the localized subspace, which, in the language of cRPA, reduces the screening from the rest of the system.  I. Calculated values of the on-site (U ) and intersite (V0i) interactions for graphite, graphene and benzene. Values are given in eV. The notation V0i denotes the intersite interaction between an atom in the unit cell, and its i-th neighbor.
In the case of graphite, two values are indicated corresponding to the two sublattices of the system.
Already in the 1950s, in the context of π-conjugated systems, it was proposed by Pariser, Parr, and Pople a one-band model with nearest-neighbor hopping and intersite interaction. This model has been widely studied and few expressions have been proposed to interpolate the Coulomb interaction between the 1/r behavior at long distance, and the short range on-site value, which is the Hubbard U . In order to get a more physical insight on the values obtained by our method, we compare them to the popular Ohno interpolation formula, which reads In this expression, the intersite interaction between atoms i and j separated by the distance r ij is estimated from the on-site interaction U and an effective dielectric constant .
The fact that we can fit the values of the intersite interaction with the Ohno potential indicates that the Coulomb integrals are properly computed. The most interesting point is that the effective dielectric constant = 2.15 that is used here matches the calculated intersite interaction values matches reasonably very well with the effective dielectric constant of graphite of 2.5 found experimentally or from cRPA calculations 21 . This shows that the functional correctly describe the screening at place in graphite. Fig. 4 shows the band structure of graphene using the extended ACBN0 functional, including both s and p orbitals, orthogonalized using the Löwdin orthonormalization, as explain above. The comparison of the band dispersion close to the K point (right panel of Fig. 4) shows that the extended ACBN0 leads to Fermi velocity quite close to the GW one around the K point, compared to the LDA calculation. This demonstrates that our functional to be the experimental one of a = 5.431Å. Overall, our results are found to be in reasonable agreement with the values of Ref. 27. However, it is worth noting that the band structure of silicon computed from LDA+U + V almost matches perfectly the one obtained from linear response values. This shows that whereas the values of U and V are not fully transferable from one implementation to another one, the observables obtained from the two approaches, the extended ACBN0 functional and linear response, give very similar results. We checked that using the generalized gradient approximation instead of the LDA one for the exchange-correlation part does not lead to a significant change of the calculated values of the effective electronic parameters.

V. CONCLUSION
In conclusion, we presented in this work an efficient method to compute ab initio and self-consistently the effective electronic parameters U , J, and V . We implemented the DFT+U + V and our novel energy functional in the real-space TDDFT code Octopus. We presented results for ground-state calculation showing that our implementation yield results in good agreement with the ones previously reported in the literature. We applied our functional to a nodal-line semimetal, ZrSiSe, showing that our functional produce very similar results obtained from a by far more expensive hybrid functional calculation. Applied to low-dimensional sp compounds, our functional gives results in qualitative agreement with cRPA calculations. Finally, we tested our functional on bulk silicon and we found that our functional reproduce well the results of linear-response in supercell 27. Let us comment on the choice of localized orbitals. In this work we employed pseudo-atomic orbitals obtained from the pseudopotentials. However, our approach is not limited to pseudopotential-based codes and can straightforwardly be used using any type of localized orbitals, such as for instance Wannier orbitals. Finally, we note that in this work we followed Ref. 27 and only considered specific intersite interaction. Determining how reliable and general is this approximation will require further investigations and would require extending the presented energy functional to include other intersite interactions. The method we presented here is, with this respect, general enough such that one could easily extend it to include other interaction terms. The extension of this functional to the time-dependent case, or to compute forces and vibrational properties of solids will be investigated in a future work.
interaction energy for noncollinear spin systems For the inter-site interaction, the double counting term is given by Putting everything together, one obtains that the rotationally-invariant form corresponding to Eq.