Through the use of the piecewise-linearity condition of the total energy, we correct the self-interaction for the study of polarons by constructing nonempirical functionals at the semilocal level of theory. We consider two functionals, the γDFT and μDFT functionals, both of which are based on the addition of a weak local potential to the semilocal Hamiltonian to enforce the piecewise-linearity condition. We show that the resulting polaron properties are in good agreement with reference hybrid functional calculations. This supports the use of semilocal functionals for calculating polaron properties.

The electron self-interaction is a spurious interaction that arises in density functional theory (DFT) due to inherent approximations in the exchange-correlation functional. For a many-electron system, two descriptions of the self-interaction can be distinguished. In an orbital-dependent formulation, the one-body self-interaction results from the self-terms of the Coulomb interaction and is found to vanish in Hartree-Fock theory, where the Hartree and exchange self-terms exactly compensate. The many-body self-interaction corresponds to the deviation of the energy functional from the piecewise-linearity condition upon electron addition.1–5 The piecewise linearity is a property of the exact functional, as first emphasized by Perdew et al.1 and by Yang et al.,4 with two complementary interpretations based on ensemble states and on replicas, respectively. Standard energy functionals deviate from the piecewise linearity, lacking derivative discontinuities at integer electron occupations.1,6 Suppressing the many-body self-interaction coincides with restoring the piecewise linearity of the total energy and thus with achieving derivative discontinuities of the energy at integer occupations.1,6

The piecewise-linearity condition holds implications for systems involving integer number of electrons. For instance, in the hydrogen cation system, a piecewise-linear functional yields equilibrium bond length and dissociation energy closer to experiment compared to standard semilocal functionals.2 Another significant aspect is the ability of a functional to localize charges, which is often not achieved by standard semilocal functionals. This is due to the electron self-interaction, which causes charge delocalization. A prototypical example of charge delocalization is represented by polarons, which are quasiparticles consisting of a localized charge coupled with self-induced lattice distortions.7 The critical role of the self-interaction in determining the charge localization of polaronic states has led to extensive studies in the literature.8–14 In semiconductors or insulators, the formation of a polaron results in an energy level within the bandgap. In particular, an electron polaron is associated with an occupied energy level below the conduction band, while a hole polaron is associated with an unoccupied energy level above the valence band. The electron self-interaction reduces the energy separation between the polaron state and the respective delocalized band edge, thus favoring charge delocalization. When the energy cost of polaronic lattice distortions overcomes the energy gain from charge localization, the polaron state destabilizes. Consequently, the polaron charge delocalizes and the polaronic distortions vanish. These aspects are crucial when modeling small polarons, which are localized over short length scales, comparable to the extent of lattice bonds, and are associated with deep energy levels in the bandgap. Suppressing the self-interaction, thus, becomes essential in polaron calculations.

Various approaches have been proposed to address the self-interaction in DFT. In 1981, Perdew and Zunger introduced a method to suppress the self-interaction for each electronic state.15 Subsequently, approaches that specifically address the self-interaction of excess charges were proposed.16,17 More recently, Sio et al. introduced an ab initio formulation for polarons and derived therefrom a self-interaction-corrected functional.11,12 However, it remains unclear to which extent these approaches guarantee the piecewise linearity of the energy functional. Meta-GGA functionals partially tackle the absence of the derivative discountinuity but may struggle with localizing polarons.18 Hybrid functionals, which admix a fraction of Fock exchange to the semilocal exchange, can be tuned to suppress the self-interaction leading to polaron localization.6,8,10,19–35 However, the computational cost of hybrid functional calculations significantly exceeds that of semilocal calculations, making them impractical for large systems. Hubbard-corrected DFT+ U functionals can suppress the many-body self-interaction by introducing on-site Coulomb repulsion for the occupation of specific orbitals.9,23,24,36–46 In the context of polarons,9,23,24,46 these orbitals are selected in consideration of their contribution to the polaron state. In another method, the piecewise-linearity condition is enforced on each electron state.47,48 Among previous approaches to address the self-interaction at the semilocal level, it is also worth mentioning a scheme relying on the electron chemical potential as slope of the total energy.27 

The γDFT semilocal functional49,50 corrects for the self-interaction error through the addition of a weak local potential solely dependent on the polaron density rather than on specific orbitals. The polaron properties obtained with hybrid functionals, DFT+ U, and γDFT have been extensively compared.46,49–51 Notably, it has been observed that enforcing the piecewise-linearity condition offers the advantage of yielding ground state and transport properties of polarons that remain robust irrespective of the adopted functional.46,49–51 Moreover, this robustness becomes pivotal for the reliability of theoretical predictions, particularly in the case of polarons, which often lack experimental data.

When addressing the self-interaction, it is important to distinguish between one-electron and many-electron systems. For a one-electron system, the exact functional is free from one-electron self-interaction.6,15 For a many-electron system, one-body and many-body forms of self-interaction can be defined. When using hybrid functionals or semilocal functionals, the description of the many-body self-interaction is in competition with the description of the one-body self-interaction.49 To illustrate this, we consider the global hybrid functional PBE0( α),52 which admixes a fraction α of Fock exchange with a fraction 1 α of Perdew–Burke–Ernzerhof (PBE) exchange.53 For α = 1, the PBE0( α) functional reduces to a Hartree–Fock-like functional, thus suppressing the one-body self-interaction. At variance, for a fraction α = α k of Fock exchange, the total energy becomes piecewise linear, thereby suppressing the many-body self-interaction.49,50 Within this context, it has been shown that for a many-electron system, one should correct for the many-body self-interaction, as this allows one to include electron screening effects, which are not accounted for when correcting for the one-body self-interaction.49,50 Specifically, for the PBE0( α) hybrid functional, the many-body self-interaction energy correction can be expressed as the one-body self-interaction energy correction divided by the high-frequency dielectric constant.49,50

In this work, we construct nonempirical semilocal functionals to model polarons by suppressing the many-body self-interaction in the polaron state. We here introduce the μDFT functional, which is conceptually similar to the γDFT functional developed previously.49,50 In both functionals, a weak local potential, which is linearly dependent on the polaron density, is added to the semilocal Hamiltonian. In the γDFT functional, the localizing potential is constructed as a derivative of the exchange-correlation functional with respect to the polaron charge, modulated by a strength γ. In the μDFT functional, the localizing potential is taken to be proportional to the polaron density modulated by a strength μ. The parameters γ and μ are nonempirically fixed by enforcing the piecewise-linearity condition of the total energy. This is achieved by imposing the condition of constant polaron level upon charge addition, after the inclusion of proper finite-size corrections in the energetics. We show that both functionals yield polaron properties in agreement with reference hybrid functional calculations, despite their semilocal functional form. This supports the use of nonempirical semilocal functionals to study polarons in DFT calculations.

We illustrate our schemes for localizing polarons at the semilocal level of theory. The fundamental concept behind these schemes is that polaron localization can be achieved by introducing a weak local potential to the semilocal Hamiltonian, promoting localization over delocalization. This stems from the observation that localized and delocalized states are often in close energetic competition.54,55

The notion that a localized potential could effectively address the self-interaction by enforcing the piecewise-linearity condition has previously been employed in the DFT+ U scheme45 and in the approach proposed by Lany and Zunger.9 In both these schemes, the U correction is applied to specific atomic orbitals. In contrast, the γDFT and μDFT functionals are based on adding to the semilocal Hamiltonian a weak local potential that self-consistently originates from the polaron density and affects all states of the system, rather than being limited to selected atomic orbitals.

In the γDFT functional,49,50 the weak local potential is essentially the first derivative of the exchange-correlation potential with respect to charge addition, modulated by a strength γ. This results in the following expression:
H σ γ = H σ 0 + q γ V x c σ q ,
(1)
where H σ 0 is the PBE Hamiltonian,53, V x σ is the PBE exchange-correlation potential, q is the polaron charge, and σ is the spin. The total energy corresponding to the Hamiltonian in Eq. (1) is51 
E γ = E 0 + q 2 2 γ σ d r V x c σ q d n σ γ ( r ) d q ,
(2)
where E 0 is the semilocal PBE energy and n σ γ is the total electron density in the spin channel σ. In Eqs. (1) and (2), the derivative of the exchange-correlation potential is calculated by taking the finite difference of two exchange-correlation potentials, which are constructed by including or excluding the polaron density in the total density at fixed valence wave functions.
Inspired by the same idea, we introduce an alternative semilocal functional, denoted μDFT, in which the added potential to the PBE Hamiltonian depends directly on the polaron density, modulated by a strength μ, namely,
H σ μ = H σ 0 + q μ n p μ δ σ , σ p ,
(3)
where n p μ is the polaron density, σ p is the spin channel containing the polaron, and δ σ , σ p is the Kronecker delta. The parameter μ carries dimensions to ensure that the added term in Eq. (3) is a potential. The polaron density is obtained by taking the square modulus of the polaron wave function, which corresponds to the highest-occupied state in the case of electron polarons and to the lowest unoccupied state in the case of hole polarons. The total energy corresponding to the Hamiltonian in Eq. (3) is
E μ = E 0 + q 2 2 μ d r n p μ ( r ) d n σ p μ ( r ) d q .
(4)

The γDFT and μDFT functionals share several properties. First, electron and hole polarons are treated analogously since the localizing potential acts on all states in the same way. Moreover, for q = 0, the localizing potential vanishes and the electronic structure coincides with that of a regular PBE calculation. In Eqs. (2) and (4), the response of the total electron density upon the addition of the polaron charge is calculated by finite differences, which requires performing two successive calculations. The γDFT and μDFT schemes also have several practical advantages. First, the same Hamiltonian is applied to all the electronic orbitals of the system, which ensures that the orthogonality of the wave functions can be achieved through the use of standard diagonalization algorithms. Second, the localized nature of the additional potential facilitates a smooth convergence of the Kohn–Sham equations. Third, it is rather straightforward to implement a local potential in existing DFT codes.

We now discuss how to enforce the piecewise-linearity condition when using γDFT and μDFT functionals. Given the similarity of these functionals, we use ξ to commonly denote the parameters γ and μ. The piecewise-linearity condition can be enforced by setting ξ = ξ k such that the second derivative of the total energy with respect to q vanishes. Through Janak’s theorem,56 this corresponds to a polaron energy level ϵ \, p ξ being constant with respect to fractional electron occupation, namely,
d 2 d q 2 E ξ ( q ) | ξ = ξ k = J a n a k d d q ϵ p ξ ( q ) | ξ = ξ k = 0.
(5)
Under the assumption of the total energy being quadratic with respect to the polaron charge,50 the condition in Eq. (5) is achieved when the neutral and charged polaron levels coincide, namely, ϵ p ξ k ( q ) = ϵ p ξ k ( 0 ). We here illustrate the workflow to find ξ k in the case of an electron polaron, which corresponds to the addition of an electron ( q = 1). The case of a hole polaron is analogous and corresponds to the removal of an electron ( q = + 1). First, one considers an initial structure in which the symmetry of the lattice is broken to accomodate the polaron lattice distortions. This is a crucial step since localized and delocalized polaronic states can be found in close energetic competition. For instance, for an electron polaron, the energy level of the extra electron coincides with the minimum of the conduction band if the symmetry of the lattice is not broken. In this case, the energy gain due to electron localization vanishes, and a structural relaxation would yield the delocalized state. The distortion can be imagined on the basis of how the lattice would respond to polaron charge localization at a given site. Once the structure is distorted, one makes an initial guess of the parameter ξ to localize the polaron. For the chosen value of the parameter ξ, one performs structural and electronic relaxations. We note that the evaluation of the energy in Eq. (2) is not required in the search of the polaronic structure. Hence, the costs of γDFT and μDFT structural relaxations essentially coincide with that of PBE calculations. For a sufficiently large value of ξ, the relaxed structure is characterized by polaronic structural distortions, which we denote R 1 ξ. Given the structure R 1 ξ, one evaluates the extent by which the piecewise-linearity condition is satisfied. The polaron level ϵ p ξ ( 1 ) is found as the Kohn–Sham level of the last occupied state in a calculation with q = 1 and structure R 1 ξ. The polaron level ϵ p ξ ( 0 ) is found as the Kohn–Sham level of the first unoccupied state in a calculation with q = 0 and structure R 1 ξ. We note that the neutral energy level ϵ p ξ ( 0 ) is equal to the PBE level ϵ p 0 ( 0 ) because of the prefactor q in Eqs. (1) and (3). Then, one repeats the same procedure for another value of ξ that leads to polaronic lattice distortions. Assuming the linearity of polaron levels with ξ,50 one can find an initial estimate of ξ k, such that ϵ p ξ k ( 1 ) = ϵ p ξ k ( 0 ). The structure is then optimized for such value of ξ k, and the procedure is repeated until ξ k converges.
In the search of ξ k, the polaron level ϵ p ξ ( q ) could be close to that of a delocalized band state. In particular, the electron polaron level could overlap with the valence band levels. Similarly, the hole polaron level could overlap with the conduction band levels. This issue can be overcome through the use of a self-consistent scissor operator,51 which is added to the Hamiltonian to shift the energy levels of the band resonating with the polaron state. Such a scissor operator can be constructed as a sum of projectors defined with the wave functions ψ i σ ξ obtained during the self-consistent optimization of the Kohn–Sham equations, namely,
S σ = Δ i M σ ξ | ψ i σ ξ ψ i σ ξ | ,
(6)
where Δ is the energy shift and M σ ξ is the manifold of band states affected by the scissor operator, namely, the valence band manifold for electron polarons and the conduction band manifold for hole polarons. The inclusion of the scissor operator S σ in the Hamiltonian shifts the total energy by a contribution N Δ, where N is the number of valence electrons.
The polaron stability can be measured through the concept of formation energy. This is defined as57,
E f ξ ( q ) = E ξ ( q ) E ref ξ ( 0 ) + q ϵ b ξ ,
(7)
where E ξ ( q ) is the total energy of the polaron system, E ref ξ ( 0 ) is the total energy of the neutral undistorted system, and ϵ b ξ is the band level corresponding to the delocalized state. For electron polarons, ϵ b ξ is the energy level of the conduction band minimum, while for hole polarons, ϵ b ξ is the energy level of the valence band maximum. When the piecewise linearity condition holds [Eq. (5)], the total energy can be rewritten as50,
E ξ k ( q ) = E 0 ( 0 ) q ϵ p ξ k ,
(8)
which leads to the following expression for the polaron formation energy:50 
E f ξ k ( q ) = q ( ϵ b ξ k ϵ p ξ k ) + [ E ξ k ( 0 ) E ref ξ k ( 0 ) ] .
(9)
For the formation energy in Eq. (9), the calculation of the total energies in Eqs. (2) and (4) is not required, making the computational cost of γDFT and μDFT calculations effectively equivalent to that of a PBE calculation. Additionally, we remark that the expression in Eq. (9) carries a transparent physical meaning. Indeed, q ( ϵ b ξ k ϵ p ξ k ) is the energy of a vertical excitation from the polaron state to the delocalized state, while E ξ k ( 0 ) E ref ξ k ( 0 ) is the cost of the involved lattice distortions.
The use of supercells for modeling polarons implies spurious finite-size effects in the energetics. Indeed, the electrostatic potential generated by a polaronic defect interacts with that of periodic replicas and with the background charge in the system, thereby affecting the total energy and the polaron energy levels.58 These effects are present for both the charged and neutral polaronic states.59 To overcome these issues, finite-size electrostatic corrections need to be applied.57–60 This can be achieved using the scheme introduced by Falletta, Wiktor, and Pasquarello,59 which extends the finite-size correction scheme of Freysoldt, Neugebauer, and Van de Walle60 to the case of defects involving built-in lattice distortions. The correction for the total energy of a calculation with polaronic charge q is59,
E cor ( q , R q ) = E m ( q , ε 0 ) E m ( q + q pol , ε ) + E m ( q + q pol , ε ) ,
(10)
where E m ( q , ε ) is the correction to the total energy of a system of charge q and screening described by a dielectric constant ε,58,60 ε is the high-frequency dielectric constant, and ε 0 is the static dielectric constant, and q pol = q ( 1 ε / ε 0 ) is the ionic polarization charge associated to the polaron lattice distortions. The charge q is set to 1 for electron polarons and to + 1 for hole polarons. The charge q is set to 0 for the neutral polaronic state, to 1 for the electron polaron state, and to + 1 for the hole polaron state. The correction for the polaron level is obtained through Janak’s theorem and is given by59,
ϵ cor ( q , R q ) = 2 E m ( q + q pol , ε ) q + q pol .
(11)
The code for performing these corrections is freely available.61,62 Without the inclusion of such finite-size corrections, the polaron formation energies would be noticeably underestimated.46,50 For simplicity, we here include these corrections in the calculated energetics, without explicitly showing them in the equations.

The calculations are performed with the quantum espresso suite.63 The implementations of the γDFT functional and of the self-consistent scissor operator are available in version 7.2, while the implementation of the μDFT functional will be made available for incorporation in the next official release of the code. The core–valence interactions are described by normconserving pseudopotentials.64 We consider the electron polaron in BiVO 4 and the hole polaron in MgO. We model BiVO 4 with a 96-atom orthorhombic supercell ( a = 10.34 Å, b = 10.34 Å, c = 11.79 Å) and MgO with a 64-atom cubic supercell ( a = 8.45 Å). The energy cutoff is set to 100 Ry in all calculations. The lattice parameters are determined at the PBE level of theory for the pristine system. The Brillouin zone is sampled at the Γ point. Convergence tests involving larger simulations cells and denser k-point samplings have shown that the polaron formation energies obtained with the present computational setup are accurate within 0.06 eV.50 Through the application of finite electric fields65 at the PBE level of theory, we determine the high-frequency and static dielectric constants, which are used for the finite-size corrections. For BiVO 4, our calculations give ε = 5.83 and ε 0 = 64.95. For MgO, we obtain ε = 2.77 and ε 0 = 10.73. In the following, we perform polaron calculations with the γDFT and μDFT semilocal functionals and compare them with reference PBE0( α) hybrid functional calculations.49,50

First, we determine the parameters of the functionals such that the piecewise-linearity condition in Eq. (5) is satisfied. With this aim, we initialize the structure with polaronic distortions that break the symmetry of the lattice. For the electron polaron in BiVO 4, we elongate the V–O bonds around a V atom. For the hole polaron in MgO, we elongate the O–Mg bonds around an O atom. Next, we perform electronic and structural relaxations to find the parameter ξ k and the corresponding structure R q ξ k, as described in Sec. II. The parameters found with this procedure are given in Table I. As an example, we show in Fig. 1 the charged and neutral polaron levels obtained with the μDFT functional at the geometry R q ξ k. We remark that the neutral level ϵ p μ ( 0 ) is essentially independent of μ since the Hamiltonian H σ μ ( 0 ) depends on μ only through the structure R q μ, which varies only marginally with μ.

FIG. 1.

Polaron energy levels ϵ p μ ( q ) and ϵ p μ ( 0 ) for the structure R q μ k as a function of μ for (a) the electron polaron in BiVO 4 and (b) the hole polaron in MgO. The polaron levels are identified by their respective polaron charge. The value μ k is found such that ϵ p μ k ( q ) = ϵ p μ k ( 0 ). The charge q is 1 for electron polarons and + 1 for hole polarons.

FIG. 1.

Polaron energy levels ϵ p μ ( q ) and ϵ p μ ( 0 ) for the structure R q μ k as a function of μ for (a) the electron polaron in BiVO 4 and (b) the hole polaron in MgO. The polaron levels are identified by their respective polaron charge. The value μ k is found such that ϵ p μ k ( q ) = ϵ p μ k ( 0 ). The charge q is 1 for electron polarons and + 1 for hole polarons.

Close modal
TABLE I.

The parameters γk, μk, and αk that enforce the piecewise-linearity condition on the polaron state for the γDFT, μDFT, and PBE0(α) functionals, respectively. The values of μk are given in Ry⋅bohr−3.

Systemγkμkαk
BiVO4 1.80 2.85 0.14 
MgO 1.96 10.90 0.34 
Systemγkμkαk
BiVO4 1.80 2.85 0.14 
MgO 1.96 10.90 0.34 
Having selected the parameter of the functional in a nonempirical way to suppress the many-body self-interaction, we now compare the electronic and structural properties of the obtained polarons. In particular, we consider the polaron densities obtained with γDFT, μDFT, and PBE0( α k). In Fig. 2, we show these densities integrated over the x y planes, namely,
n p ( z ) = d x d y n p ( x , y , z )
(12)
and find an excellent agreement among the three considered piecewise-linear functionals. We then consider the polaronic lattice distortions and find that the polaronic bonds obtained with γDFT and μDFT essentially concide and are in very good agreement with PBE0( α k) results, with deviations of at most 0.02 Å (cf. Table II).
FIG. 2.

Polaron densities for (a) the electron polaron in BiVO 4 and (b) the hole polaron in MgO as obtained with the piecewise-linear functionals γDFT, μDFT, and PBE0( α k). The polaron densities are integrated over x y-planes.

FIG. 2.

Polaron densities for (a) the electron polaron in BiVO 4 and (b) the hole polaron in MgO as obtained with the piecewise-linear functionals γDFT, μDFT, and PBE0( α k). The polaron densities are integrated over x y-planes.

Close modal
TABLE II.

V–O and Mg–O bond lengths (in Å) of the distorted polaronic structures as obtained with the piecewise-linear functionals γDFT, μDFT, and PBE0(αk) for electron and hole polarons in BiVO4 and MgO, respectively.

SystemγDFTμDFTPBE0 (αk)
BiVO4 1.82 1.82 1.80 
MgO 2.23 2.23 2.20 
SystemγDFTμDFTPBE0 (αk)
BiVO4 1.82 1.82 1.80 
MgO 2.23 2.23 2.20 
We verify that the localizing potential in μDFT is indeed weak. For this purpose, we compare this potential with the Hartree potential generated by the polaron charge density. In Fig. 3, we show these potentials averaged over x y planes as a function of the z coordinate, i.e.,
V ( z ) = 1 A x y d x d y V ( x , y , z ) ,
(13)
where V is the potential, and A x y is the area of x y planes in the supercell. We see that the added localizing potential in μDFT is weak and carries an opposite sign with respect to the electrostatic potential, thus favoring polaron localization. The fact that the added potential in μDFT is more localized than the Hartree potential generated by the polaron density facilitates the self-consistent solution of the Kohn–Sham equations.
FIG. 3.

Electrostatic potential V elec = V H [ q n p μ k ] and μDFT localizing potential V σ p μ k = q μ k n p μ k averaged over x y-planes for (a) the electron polaron in BiVO 4 and (b) the hole polaron in MgO. The charge q is 1 for electron polarons and + 1 for hole polarons.

FIG. 3.

Electrostatic potential V elec = V H [ q n p μ k ] and μDFT localizing potential V σ p μ k = q μ k n p μ k averaged over x y-planes for (a) the electron polaron in BiVO 4 and (b) the hole polaron in MgO. The charge q is 1 for electron polarons and + 1 for hole polarons.

Close modal

Finally, we compare the formation energies obtained with the piecewise-linear functionals using the expression in Eq. (9). We find that the formation energies obtained with γDFT and μDFT are very close, with deviations of at most 0.03 eV (cf. Table III). A good agreement is also found with the PBE0( α k) results. This further corroborates the finding that upon enforcing the piecewise-linearity condition the polaron properties turn out to weakly depend on the adopted functional.46,49–51 This is particularly relevant in the case of the hole polaron in MgO, for which the experimental characterization remains inconclusive66,67 and the theoretical description contentious.8,10,50,59,68

TABLE III.

Polaron formation energies (in eV) obtained with the piecewise-linear functionals γDFT, μDFT, and PBE0(αk). A negative formation energy indicates polaron stability.

γDFTμDFTPBE0(αk)
BiVO4 −0.44 −0.44 −0.63 
MgO −0.50 −0.47 −0.53 
γDFTμDFTPBE0(αk)
BiVO4 −0.44 −0.44 −0.63 
MgO −0.50 −0.47 −0.53 

We investigate in this work the performance of suitably constructed semilocal density functionals in addressing polaron properties. As case studies, we consider both electron and hole polarons, for which our approach works in an equivalent manner. The two considered semilocal functionals, denoted γDFT and μDFT, depend on a single parameter, which we fix in a nonempirical fashion by enforcing the piecewise-linearity condition and thus suppressing the many-body self-interaction. Both functionals are successful in achieving polaron localization. We find that the achieved results are in excellent agreement with those of reference hybrid-functional calculations. The present analysis encompasses both electronic and structural polaron properties, such as polaron densities, polaronic bond distortions, and polaron formation energies. This corroborates the finding that these polaron properties are robust and that they only weakly depend on the functional adopted, provided the piecewise-linearity condition is satisfied. Overall, our study shows that polaron properties can be accurately described through the use of semilocal functionals, without resorting to computationally more expensive hybrid-functional calculations.

We thank Michele Pavanello for useful discussions. The calculations have been performed at the Swiss National Supercomputing Centre (CSCS) (Grant under Projects ID No. s1122) and at SCITAS-EPFL.

The authors declare no conflict of interest.

Stefano Falletta: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Alfredo Pasquarello: Conceptualization (equal); Data curation (supporting); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (lead); Supervision (lead); Validation (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

The data that support the findings of this study are openly available in Materials Cloud, Ref. 69.

1.
J. P.
Perdew
,
R. G.
Parr
,
M.
Levy
, and
J. L.
Balduz
, “
Density-functional theory for fractional particle number: Derivative discontinuities of the energy
,”
Phys. Rev. Lett.
49
,
1691
(
1982
).
2.
A.
Ruzsinszky
,
J. P.
Perdew
,
G. I.
Csonka
,
O. A.
Vydrov
, and
G. E.
Scuseria
, “
Density functionals that are one- and two- are not always many-electron self-interaction-free, as shown for H 2 +, He 2 +, LiH +, and Ne 2 +
,”
J. Chem. Phys.
126
,
104102
(
2007
).
3.
Y.
Zhang
and
W.
Yang
, “
A challenge for density functionals: Self-interaction error increases for systems with a noninteger number of electrons
,”
J. Chem. Phys.
109
,
2604
(
1998
).
4.
W.
Yang
,
Y.
Zhang
, and
P. W.
Ayers
, “
Degenerate ground states and a fractional number of electrons in density and reduced density matrix functional theory
,”
Phys. Rev. Lett.
84
,
5172
(
2000
).
5.
P.
Mori-Sánchez
,
A. J.
Cohen
, and
W.
Yang
, “
Many-electron self-interaction error in approximate density functionals
,”
J. Chem. Phys.
125
,
201102
(
2006
).
6.
L.
Kronik
and
S.
Kümmel
, “
Piecewise linearity, freedom from self-interaction, and a Coulomb asymptotic potential: Three related yet inequivalent properties of the exact density functional
,”
Phys. Chem. Chem. Phys.
22
,
16467
(
2020
).
7.
C.
Franchini
,
M.
Reticcioli
,
M.
Setvin
, and
U.
Diebold
, “
Polarons in materials
,”
Nat. Rev. Mater
6
,
560
(
2021
).
8.
J. B.
Varley
,
A.
Janotti
,
C.
Franchini
, and
C. G.
Van de Walle
, “
Role of self-trapping in luminescence and p-type conductivity of wide-band-gap oxides
,”
Phys. Rev. B
85
,
081109
(
2012
).
9.
S.
Lany
and
A.
Zunger
, “
Polaronic hole localization and multiple hole binding of acceptors in oxide wide-gap semiconductors
,”
Phys. Rev. B
80
,
085202
(
2009
).
10.
S.
Kokott
,
S. V.
Levchenko
,
P.
Rinke
, and
M.
Scheffler
, “
First-principles supercell calculations of small polarons with proper account for long-range polarization effects
,”
New J. Phys.
20
,
033023
(
2018
).
11.
W. H.
Sio
,
C.
Verdi
,
S.
Poncé
, and
F.
Giustino
, “
Ab initio theory of polarons: Formalism and applications
,”
Phys. Rev. B
99
,
235139
(
2019
).
12.
W. H.
Sio
,
C.
Verdi
,
S.
Poncé
, and
F.
Giustino
, “
Polarons from first principles, without supercells
,”
Phys. Rev. Lett.
122
,
246403
(
2019
).
13.
J.
Lafuente-Bartolome
,
C.
Lian
,
W. H.
Sio
,
I. G.
Gurtubay
,
A.
Eiguren
, and
F.
Giustino
, “
Unified approach to polarons and phonon-induced band structure renormalization
,”
Phys. Rev. Lett.
129
,
076402
(
2022
).
14.
V.
Vasilchenko
,
A.
Zhugayevych
, and
X.
Gonze
, “
Variational polaron equations applied to the anisotropic Fröhlich model
,”
Phys. Rev. B
105
,
214301
(
2022
).
15.
J. P.
Perdew
and
A.
Zunger
, “
Self-interaction correction to density-functional approximations for many-electron systems
,”
Phys. Rev. B
23
,
5048
(
1981
).
16.
M.
d’Avezac
,
M.
Calandra
, and
F.
Mauri
, “
Density functional theory description of hole-trapping in SiO 2: A self-interaction-corrected approach
,”
Phys. Rev. B
71
,
205210
(
2005
).
17.
J.
VandeVondele
and
M.
Sprik
, “
A molecular dynamics study of the hydroxyl radical in solution applying self-interaction-corrected density functional methods
,”
Phys. Chem. Chem. Phys.
7
,
1363
(
2005
).
18.
D.
Wickramaratne
and
J. L.
Lyons
, “Shortcomings of using the scan functional for point defects and polarons in semiconductors,” arXiv:2311.03634 [cond-mat.mtrl-sci] (2023).
19.
A. R.
Elmaslmane
,
J.
Wetherell
,
M. J. P.
Hodgson
,
K. P.
McKenna
, and
R. W.
Godby
, “
Accuracy of electron densities obtained via Koopmans-compliant hybrid functionals
,”
Phys. Rev. Mater.
2
,
040801
(
2018
).
20.
J. J.
Carey
and
K. P.
McKenna
, “
Screening doping strategies to mitigate electron trapping at anatase TiO 2 surfaces
,”
J. Phys. Chem. C
123
,
22358
(
2019
).
21.
J. A.
Quirk
,
V. K.
Lazarov
, and
K. P.
McKenna
, “
First-principles modeling of oxygen-deficient anatase TiO 2 nanoparticles
,”
J. Phys. Chem. C
124
,
23637
(
2020
).
22.
J. J.
Carey
,
J. A.
Quirk
, and
K. P.
McKenna
, “
Hole polaron migration in bulk phases of TiO 2 using hybrid density functional theory
,”
J. Phys. Chem. C
125
,
12441
(
2021
).
23.
J. J.
Carey
and
K. P.
McKenna
, “
Does polaronic self-trapping occur at anatase TiO 2 surfaces?
,”
J. Phys. Chem. C
122
,
27540
(
2018
).
24.
A. R.
Elmaslmane
,
M. B.
Watkins
, and
K. P.
McKenna
, “
First-principles modeling of polaron formation in TiO 2 polymorphs
,”
J. Chem. Theory Comput.
14
,
3740
(
2018
).
25.
G.
Miceli
,
W.
Chen
,
I.
Reshetnyak
, and
A.
Pasquarello
, “
Nonempirical hybrid functionals for band gaps and polaronic distortions in solids
,”
Phys. Rev. B
97
,
121112
(
2018
).
26.
P.
Deák
,
Q.
Duy Ho
,
F.
Seemann
,
B.
Aradi
,
M.
Lorke
, and
T.
Frauenheim
, “
Choosing the correct hybrid for defect calculations: A case study on intrinsic carrier trapping in β-Ga 2O 3
,”
Phys. Rev. B
95
,
075208
(
2017
).
27.
B.
Sadigh
,
P.
Erhart
, and
D.
Åberg
, “
Variational polaron self-interaction-corrected total-energy functional for charge excitations in insulators
,”
Phys. Rev. B
92
,
075202
(
2015
).
28.
N.
Sai
,
P. F.
Barbara
, and
K.
Leung
, “
Hole localization in molecular crystals from hybrid density functional theory
,”
Phys. Rev. Lett.
106
,
226403
(
2011
).
29.
S.
Refaely-Abramson
,
S.
Sharifzadeh
,
M.
Jain
,
R.
Baer
,
J. B.
Neaton
, and
L.
Kronik
, “
Gap renormalization of molecular crystals from density-functional theory
,”
Phys. Rev. B
88
,
081204
(
2013
).
30.
T.
Bischoff
,
J.
Wiktor
,
W.
Chen
, and
A.
Pasquarello
, “
Nonempirical hybrid functionals for band gaps of inorganic metal-halide perovskites
,”
Phys. Rev. Mater.
3
,
123802
(
2019
).
31.
T.
Bischoff
,
I.
Reshetnyak
, and
A.
Pasquarello
, “
Adjustable potential probes for band-gap predictions of extended systems through nonempirical hybrid functionals
,”
Phys. Rev. B
99
,
201114
(
2019
).
32.
N.
Österbacka
,
P.
Erhart
,
S.
Falletta
,
A.
Pasquarello
, and
J.
Wiktor
, “
Small electron polarons in CsPbBr 3: Competition between electron localization and delocalization
,”
Chem. Mater.
32
,
8393
(
2020
).
33.
T.
Bischoff
,
I.
Reshetnyak
, and
A.
Pasquarello
, “
Band gaps of liquid water and hexagonal ice through advanced electronic-structure calculations
,”
Phys. Rev. Res.
3
,
023182
(
2021
).
34.
J.
Yang
,
S.
Falletta
, and
A.
Pasquarello
, “
One-shot approach for enforcing piecewise linearity on hybrid functionals: Application to band gap predictions
,”
J. Phys. Chem. Lett.
13
,
3066
(
2022
).
35.
J.
Yang
,
S.
Falletta
, and
A.
Pasquarello
, “
Range-separated hybrid functionals for accurate prediction of band gaps of extended systems
,”
npj Comput. Mater.
9
,
108
(
2023
).
36.
V. I.
Anisimov
and
O.
Gunnarsson
, “
Density-functional calculation of effective Coulomb interactions in metals
,”
Phys. Rev. B
43
,
7570
(
1991
).
37.
V. I.
Anisimov
,
J.
Zaanen
, and
O. K.
Andersen
, “
Band theory and Mott insulators: Hubbard U instead of stoner I
,”
Phys. Rev. B
44
,
943
(
1991
).
38.
V. I.
Anisimov
,
I. V.
Solovyev
,
M. A.
Korotin
,
M. T.
Czyżyk
, and
G. A.
Sawatzky
, “
Density-functional theory and NiO photoemission spectra
,”
Phys. Rev. B
48
,
16929
(
1993
).
39.
I. V.
Solovyev
,
P. H.
Dederichs
, and
V. I.
Anisimov
, “
Corrected atomic limit in the local-density approximation and the electronic structure of d impurities in Rb
,”
Phys. Rev. B
50
,
16861
(
1994
).
40.
M. T.
Czyżyk
and
G. A.
Sawatzky
, “
Local-density functional and on-site correlations: The electronic structure of La 2CuO 4 and LaCuO 3
,”
Phys. Rev. B
49
,
14211
(
1994
).
41.
A. I.
Liechtenstein
,
V. I.
Anisimov
, and
J.
Zaanen
, “
Density-functional theory and strong interactions: Orbital ordering in Mott-Hubbard insulators
,”
Phys. Rev. B
52
,
R5467
(
1995
).
42.
V. I.
Anisimov
,
F.
Aryasetiawan
, and
A. I.
Lichtenstein
, “
First-principles calculations of the electronic structure and spectra of strongly correlated systems: The LDA+ U method
,”
J. Phys.: Condens. Matter
9
,
767
(
1997
).
43.
S. L.
Dudarev
,
G. A.
Botton
,
S. Y.
Savrasov
,
C. J.
Humphreys
, and
A. P.
Sutton
, “
Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+ U study
,”
Phys. Rev. B
57
,
1505
(
1998
).
44.
A. G.
Petukhov
,
I. I.
Mazin
,
L.
Chioncel
, and
A. I.
Lichtenstein
, “
Correlated metals and the LDA+ U method
,”
Phys. Rev. B
67
,
153106
(
2003
).
45.
M.
Cococcioni
and
S.
de Gironcoli
, “
Linear response approach to the calculation of the effective interaction parameters in the LDA+ U method
,”
Phys. Rev. B
71
,
035105
(
2005
).
46.
S.
Falletta
and
A.
Pasquarello
, “
Hubbard U through polaronic defect states
,”
npj Comput. Mater.
8
,
265
(
2022
).
47.
I.
Dabo
,
A.
Ferretti
,
N.
Poilvert
,
Y.
Li
,
N.
Marzari
, and
M.
Cococcioni
, “
Koopmans’ condition for density-functional theory
,”
Phys. Rev. B
82
,
115121
(
2010
).
48.
N. L.
Nguyen
,
N.
Colonna
,
A.
Ferretti
, and
N.
Marzari
, “
Koopmans-compliant spectral functionals for extended systems
,”
Phys. Rev. X
8
,
021051
(
2018
).
49.
S.
Falletta
and
A.
Pasquarello
, “
Many-body self-interaction and polarons
,”
Phys. Rev. Lett.
129
,
126401
(
2022
).
50.
S.
Falletta
and
A.
Pasquarello
, “
Polarons free from many-body self-interaction in density functional theory
,”
Phys. Rev. B
106
,
125119
(
2022
).
51.
S.
Falletta
and
A.
Pasquarello
, “
Polaron hopping through piecewise-linear functionals
,”
Phys. Rev. B
107
,
205125
(
2023
).
52.
J. P.
Perdew
,
M.
Ernzerhof
, and
K.
Burke
, “
Rationale for mixing exact exchange with density functional approximations
,”
J. Chem. Phys.
105
,
9982
(
1996
).
53.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
(
1996
).
54.
J. L.
Gavartin
,
P. V.
Sushko
, and
A. L.
Shluger
, “
Modeling charge self-trapping in wide-gap dielectrics: Localization problem in local density functionals
,”
Phys. Rev. B
67
,
035108
(
2003
).
55.
A.
Carvalho
,
A.
Alkauskas
,
A.
Pasquarello
,
A. K.
Tagantsev
, and
N.
Setter
, “
A hybrid density functional study of lithium in ZnO: Stability, ionization levels, and diffusion
,”
Phys. Rev. B
80
,
195205
(
2009
).
56.
J. F.
Janak
, “
Proof that e n i = ϵ in density-functional theory
,”
Phys. Rev. B
18
,
7165
(
1978
).
57.
C.
Freysoldt
,
B.
Grabowski
,
T.
Hickel
,
J.
Neugebauer
,
G.
Kresse
,
A.
Janotti
, and
C. G.
Van de Walle
, “
First-principles calculations for point defects in solids
,”
Rev. Mod. Phys.
86
,
253
(
2014
).
58.
H.-P.
Komsa
,
T. T.
Rantala
, and
A.
Pasquarello
, “
Finite-size supercell correction schemes for charged defect calculations
,”
Phys. Rev. B
86
,
045112
(
2012
).
59.
S.
Falletta
,
J.
Wiktor
, and
A.
Pasquarello
, “
Finite-size corrections of defect energy levels involving ionic polarization
,”
Phys. Rev. B
102
,
041115
(
2020
).
60.
C.
Freysoldt
,
J.
Neugebauer
, and
C. G.
Van de Walle
, “
Fully Ab Initio finite-size corrections for charged-defect supercell calculations
,”
Phys. Rev. Lett.
102
,
016402
(
2009
).
61.
S.
Falletta
,
J.
Wiktor
, and
A.
Pasquarello
, see https://github.com/falletta/finite-size-corrections-defect-levels for the code to apply finite-size corrections to polaronic states.
62.
S.
Falletta
,
J.
Wiktor
, and
A.
Pasquarello
, “
Finite-size corrections of defect energy levels involving ionic polarization
,”
Materials Cloud Archive 2020
.
70
, 1 (
2020
).
63.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
, and
A.
Dal Corso
, “
QUANTUM ESPRESSO: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
64.
M.
van Setten
,
M.
Giantomassi
,
E.
Bousquet
,
M.
Verstraete
,
D.
Hamann
,
X.
Gonze
, and
G.-M.
Rignanese
, “
The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table
,”
Comput. Phys. Commun.
226
,
39
(
2018
).
65.
P.
Umari
and
A.
Pasquarello
, “
Ab initio molecular dynamics in a finite homogeneous electric field
,”
Phys. Rev. Lett.
89
,
157602
(
2002
).
66.
Z. A.
Rachko
and
J. A.
Valbis
, “
Luminescence of free and relaxed excitons in MgO
,”
Phys. Status Solidi B
93
,
161
(
1979
).
67.
S.
Dolgov
,
T.
Kärner
,
A.
Lushchik
,
A.
Maaroos
,
S.
Nakonechnyi
, and
E.
Shablonin
, “
Trapped-hole centers in MgO single crystals
,”
Phys. Solid State
53
,
1244
(
2011
).
68.
A. L.
Shluger
,
E. N.
Heifets
,
J. D.
Gale
, and
C. R. A.
Catlow
, “
Theoretical simulation of localized holes in MgO
,”
J. Phys.: Condens. Matter
4
,
5711
(
1992
).
69.
S.
Falletta
and
A.
Pasquarello
, “
Nonempirical semilocal density functionals for correcting the self-interaction of polaronic states
,”
Materials Cloud Archive 2024.
47
, 1 (
2024
).