Embedding for bulk systems using localized atomic orbitals

TL;DR

This paper introduces an orbital rotation embedding method for bulk semiconductor and insulator systems, implemented in VASP, effectively addressing defect and polaron challenges in silicon and titania.

Contribution

The paper presents a novel embedding approach based on orbital rotations, specifically designed for bulk systems, with implementation in VASP for improved defect and polaron modeling.

Findings

Effective modeling of defect structures in silicon

Successful simulation of polaron formation in titania

Implementation in VASP enhances practical applicability

Abstract

We present an embedding approach for semiconductors and insulators based on or- bital rotations in the space of occupied Kohn-Sham orbitals. We have implemented our approach in the popular VASP software package. We demonstrate its power for defect structures in silicon and polaron formation in titania, two challenging cases for conventional Kohn-Sham density functional theory.

Tables

Table 1. Table 1: Comparison of defect formation energies for different defect types in silicon. All energies in eV, errors in brackets [meV] are deviation from full hybrid benchmark calculation. DFT and embedding calculation energies represent the “single-shot” evaluation of the full hybrid energy functional using the corresponding non-hybrid PBE or embedded orbitals as indicated by an asterisk.

energies [eV] (error [meV])
Defect	Hybrid	DFT∗		Embedding∗
H	3.00	2.99	(-17)	3.01	(10)
X	3.01	3.01	(1)	3.04	(27)
C3V	3.05	3.03	(-19)	3.06	(18)
T	3.77	3.34	(-423)	3.75	(-20)
VJT	4.14	4.52	(377)	4.19	(44)
V	4.23	4.83	(599)	4.26	(25)

Table 2. Table 2: Same comparison as in Tab.  1 for an eight times larger supercell containing 512 atoms. Rightmost column shows the range of available experimental data taken from [ Exp1, ; Exp2, ; Exp3, ; Exp4, ; Exp5, ; Exp6, ; Exp7, ; Exp8, ; Exp9, ; Exp10, ].

energies [eV] (error [meV])
Defect	Hybrid	DFT∗		Embedding∗		Experiment
T	4.97	5.17	(198)	5.13	(165)
H	4.22	4.19	(38)	4.25	(28)	4.2-4.7
V	5.06	5.55	(482)	5.08	(22)	2.1-4.0

Table 3. Table 3: Energies of the distorted E dist subscript 𝐸 dist E_{\mathrm{dist}} and ideal E ideal subscript 𝐸 ideal E_{\mathrm{ideal}} lattice structure of charged rutile titania. The energy Δ ​ E Δ 𝐸 \Delta E is gained by forming a small polaron. Energies after seven iterations are given for different sizes of the embedded region. Different methods are labeled as follows. Hybrid: full hybrid functional calculation used as benchmark; DFT: direct evaluation of the energies using the PBE functional; DFT ∗ : evaluation of the hybrid energy functional using the orbitals from the PBE calculation; embedding: embedding calculations as function of the number of freeze-and-thaw cycles at fixed number of localized atomic orbitals N A subscript 𝑁 𝐴 N_{A} . After the unitary rotation to localize orbitals around the defect [Eq. ( 5 )], we optimize a subset N emb A ≤ N A superscript subscript 𝑁 emb 𝐴 subscript 𝑁 𝐴 N_{\mathrm{emb}}^{A}\leq N_{A} of orbitals in the embedded region A 𝐴 A as noted.

Method		cycles	$E_{\mathrm{dist}}$ [eV]	$E_{\mathrm{ideal}}$ [eV]	$\Delta E$ [meV]
Hybrid		-	-972.85	-972.33	514
DFT		-	-688.17	-687.81	-355
DFT∗		-	-962.81	-962.54	270
Embedding∗		1	-963.70	-963.52	174
$N_A=30$	$N_{\mathrm{emb}}^A=14$	2	-963.36	-963.12	240
		3	-963.47	-963.16	313
		4	-963.41	-962.98	426
		5	-963.45	-962.99	459
		6	-963.44	-962.98	462
		7	-963.45	-962.98	475
$N_A=30$	$N_{\mathrm{emb}}^A=6$	7	-963.40	-962.94	460
	$N_{\mathrm{emb}}^A=8$	7	-963.41	-962.95	459
	$N_{\mathrm{emb}}^A=12$	7	-963.45	-963.00	454
	$N_{\mathrm{emb}}^A=20$	7	-963.50	-963.02	474
	$N_{\mathrm{emb}}^A=30$	7	-963.59	-963.08	514
$N_A=90$	$N_{\mathrm{emb}}^A=30$	7	-963.55	-963.08	462
$N_A=6$	$N_{\mathrm{emb}}^A=6$	7	-963.38	-962.96	420

Equations

$$H|{\psi }_i⟩={\overline{\epsilon }}_i|{\psi }_i⟩,$$

$$|{\psi }_i⟩,\gamma (\vec{r},\vec{r}{}^{\prime })=\sum _{i=1}^{N_{\mathrm{tot}}}{f}_i⟨{\psi }_i|\vec{r}⟩⟨\vec{r}{}^{\prime }|{\psi }_i⟩,$$

$$W_{ki}:=⟨{\alpha }_k|{\psi }_i⟩,k=1\dots N_A,i=1\dots N.$$

$$W=U\cdot D\cdot V^{†}$$

$$V\left|\psi_{i}\right>=\left|\phi_{i}\right>,\quad\left|\left<\alpha_{j}\vphantom{\phi_{i}}\right|\left.\!\vphantom{\alpha_{j}}\phi_{i}\right>\right|\left\{\begin{array}[]{cc}\leq 1,&i\leq N_{A}\\ =0,&i>N_{A}\end{array}\right.$$

$${\epsilon }_i:=⟨{\varphi }_i|H|{\varphi }_i⟩,$$

$$\sum _{i=1}^N{\epsilon }_i=\sum _{i=1}^N{\overline{\epsilon }}_i$$

$$E=\sum _{i=1}^{N_{\mathrm{tot}}}{f}_i⟨{\varphi }_i|T|{\varphi }_i⟩+\frac{1}{2}{E}_H[\rho ]+E_{\mathrm{XC}}[{\gamma }_A,{\gamma }_B]$$

$$E_H[\rho ]=\frac{1}{2}\sum _{i,j=1}^{N_{\mathrm{tot}}}{f}_i{f}_j⟨ij|\mathcal{C}|ij⟩,$$

$$E_{\mathrm{XC}}^{AB}[{\gamma }_A,{\gamma }_B]=E_{XC}^A[{\gamma }_A]+E_{\mathrm{XC}}^B[{\gamma }_B]+E_{\mathrm{XC}}^{\mathrm{int}}[{\gamma }_A,{\gamma }_B].$$

$$E_{\mathrm{XC}}^{\mathrm{int}}[{\gamma }_A,{\gamma }_B]=E_{\mathrm{XC}}^B[{\gamma }_{\mathrm{tot}}]-E_{\mathrm{XC}}^B[{\gamma }_A]-E_{\mathrm{XC}}^B[{\gamma }_B]$$

$$E_{\mathrm{XC}}^{\mathrm{int}}[{\gamma }_A,{\gamma }_B]$$

Full text

I Introduction

Ab-initio electronic structure theory for bulk materials has experienced tremendous advances in many areas such as density functional theory Kohn ; BurkePerspective ; ChallengesDFT ; BurkeReview , improved post-DFT Nozier ; Langreth ; Miyake ; Fuchs ; Furche ; SOSEX ; RPA and, e.g., van der Waals functionals,tkatchenko12 as well as highly accurate quantum chemical MP1 ; MP2 ; PCC and Monte-Carlo approachesAlavi . However, many problems are still out of reach of an advanced theoretical description due to their size: the accurate description of, for example, defect structures requires both a highly accurate treatment of the local defect region, as well as the treatment of a large number of atoms of the environmentlibischprl . It is often challenging for a single method to meet both requirements. Embedding is therefore a suitable strategy to evercome this hurdle. Its underlying idea is to treat the local structure or, more generally, the subsystem of interest by a high-level method while treating the environment with the help of a numerically less demanding lower level method. Consistently combining different electronic structure methods within the same calculation is both the advantage and the challenge of the embedding approach NeugebauerSelfc .

Several embedding schemes have been proposed Nafziger ; Goodpaster1 ; Goodpaster2 ; Goodpaster3 ; Wesolovski1 ; Jacob ; libischreview ; Carter , relying either on some form of a local embedding potential $V_{\mathrm{emb}}(\mathbf{r})$ that mediates the interaction between the subsystem referred in the following as the cluster and the surrounding environment. More elaborate operator-based approaches Knizia1 ; Manby introduce a nonlocal embedding $\hat{V}_{\mathrm{emb}}(\mathbf{r},\mathbf{r^{\prime}})$. Typically, subsystems are treated in the presence of $V_{\mathrm{emb}}(\mathbf{r})$ [or $\hat{V}_{\mathrm{emb}}(\mathbf{r},\mathbf{r^{\prime}})$] using a high-level method, while the entire system is handled by density functional theory (DFT). The individual subsystem densities are then added to obtain an approximation for the total density of the entire system. While conceptually simpler, local embedding potentials feature the distinct disadvantage that no set of mutually orthogonal orbitals of the entire system exists. Consequently, evaluation of the total energy becomes challenging: in particular the kinetic energy needs to be approximated. Huang et al. Carter used an optimized effective potential method to recover the kinetic energy given a total electron density. Conversely, Fornace et al. presented an embedded mean-field theory ManbyMeanfield partitioning the one-particle density matrix of the system based on its basis functions. A single Hamiltonian then describes the entire system, avoiding any issues with evaluating the kinetic energy for cluster and environment separately. Additionally, this approach, by construction, allows for direct charge exchange between the cluster of interest and the environment. However, a direct extension to plane-wave basis sets used in periodic solid state computations seems challenging. Goodpaster et al. Goodpaster1 ; Goodpaster2 ; Goodpaster3 have presented a scheme relying on projection operators to ensure mutual orthogonality of orbitals belonging to different subsystems. In the present article, we present an alternative strategy to generate and maintain mutually orthogonal orbitals for the subsystems throughout the calculation. We determine Wannier-like orbitals localized within the cluster by performing unitary rotations within the subspace of fully occupied Kohn-Sham orbitals while the orthogonal complement of remaining orbitals resides within the environment Knizia1 . During the optimization cycle for the cluster involving an advanced functional, the environment orbitals remain frozen and thus orthogonality is preserved. This approach avoids the inaccuracies associated with approximating the kinetic energy.

In the present paper we demonstrate the power of our embedding scheme in a proof-of-principle calculation adressing two problems for which standard Kohn-Sham DFT is known to be inadequate: defects in silicon and polarons in titania. We use the following hierarchy of methods: the cluster is treated by the (expensive) hybrid functional PBEh while the environment is treated only by the PBE functional. We show that this embedding scheme implemented in the Vienna Ab Initio Simulation Package (VASP) is robust and efficient. We emphasize that the present embedding scheme is not limited to hybrid-DFT in DFT embeddings. Future extensions will adress the treatment of the cluster by RPA or quantum chemistry approaches.

II Technique

We partition a system into two parts: a cluster of interest $A$ with atomic sites $r_{j,A}$, ($j=1,\ldots,M_{A}$) with $M_{A}$ the number of atomic sites included in the cluster, and the surounding environment $B$, containing $M_{B}$ atomic sites $r_{j,B}$, ($j=1,\ldots,M_{B}$). In a first step, the entire system ($A+B$) is solved using a single, comparatively cheap exchange-correlation functional, e.g., PBE PBE ,

[TABLE]

yielding Kohn-Sham orbitals $\left|\psi_{i}\right>$ with orbital energies $\overline{\varepsilon}_{i}$ and the density matrix

[TABLE]

with occupation numbers $f_{i}\in[0,1]$, where the index $i=1,\ldots,N_{\mathrm{tot}}$ goes over all orbitals and physical spin. Note that we we have not included $k$-point sampling in the present ansatz, since it is not straightforward to treat the transformations at different $k$-points independently. We aim to find a unitary rotation within the subspace of fully occupied orbitals ($f_{i}=1$, $i=1,\ldots,N$) that yields a set of orbitals aligned with the atomic orbitals $\left|\alpha_{k}\right>$ localized around the atomic sites of the cluster. The index $k=1,\ldots,N_{A}$ of the atomic orbitals includes both the site index as well as radial and angular momentum quantum numbers. To this end we apply to the orbital overlap matrix W,

[TABLE]

a singular value decomposition according to

[TABLE]

with $D=\mathrm{diag}(\sigma_{i})$. The unitary matrix $V$ represents the rotation in the space of the $N$ occupied orbitals that optimally aligns $N_{A}$ of these orbitals with the atomic orbitals $\left|\alpha_{i}\right>$ keeping the remaining $N-N_{A}$ orbitals orthogonal. The singular values $\sigma_{i}$ provide a measure for the degree of overlap between the orbitals $\left|\alpha_{i}\right>$ and the rotated orbitals $\left|\phi_{i}\right>$,

[TABLE]

Orbitals with indices $i>N$ outside the space of occupied orbitals are unaffected by the rotation, $\left|\phi_{i}\right>=\left|\psi_{i}\right>,\,\forall i:N<i<N_{\mathrm{tot}}$. Using the rotated $\left|\phi_{i}\right>$, we can thus partition the occupied space into orbitals which have an overlap with the $\left|\alpha_{i}\right>$, and those who do not. Ideally, if the Kohn Sham orbitals are well covered by the atomic wavefunctions, we expect the singular values $\sigma_{i}$ to be close to 1.

After the orbital rotations, a subset $\left|\phi_{i}^{A}\right>$ of those $i=1,\ldots,N_{\mathrm{emb}}^{A}$ with $N_{\mathrm{emb}}^{A}\leq N_{A}$ can now be optimized based on a more expensive exchange-correlation (XC) functional $E_{\mathrm{XC}}^{A}$, e.g., a hybrid functional hybrid1 ; hybrid2 , while freezing the orthogonal complement of $N_{\mathrm{emb}}^{B}=N-N_{\mathrm{emb}}^{A}$ environment orbitals $\left|\phi_{i}^{B}\right>$. In general, the number of orbitals $N_{\mathrm{emb}}^{A}$ used in the embedding procedure may be smaller than the number of atomic basis functions $N_{A}$: in principle, one may choose any subset of the localized orbitals $0<N_{\mathrm{emb}}^{A}\leq N_{A}$. In practice, we sort the rotated orbitals by their singular values $\sigma_{i}$, and choose the $N_{\mathrm{emb}}^{A}$ orbitals corresponding to the largest $\sigma_{i}$, where $N_{\mathrm{emb}}^{A}$ is chosen according to the number of orbitals of interest within the cluster. A typical cut-off will be $\sigma_{i}>0.5$. We find that our results do not strongly depend on $N_{\mathrm{emb}}^{A}$, as long as the number of optimized orbitals is sufficiently large as to properly describe the local bonding. We will discuss the choice of $N_{\mathrm{emb}}^{A}$ in more detail in the results section below.

Note that after the orbital rotation, the $\left|\phi_{i}\right>$ are no longer eigenvectors of the Kohn-Sham Hamiltonian $H$. The diagonals of the Hamiltonian are given by the expectation values

[TABLE]

related to the original eigenvalues through the invariance of the trace

[TABLE]

under unitary rotations. In the case of fractional occupations of some of the orbitals in $N_{\mathrm{tot}}$, the above considerations remain valid in the subspace of the $N\leq N_{\mathrm{tot}}$ fully occupied orbitals. All orbitals with fractional occupation are assigned to cluster $A$, even if they are not well localized. However, extension of the present approach to quantum-chemistry based correlated wavefunction approaches may face difficulties for such delocalized metallic states.

Taking into consideration that different exchange-correlation functionals will be employed for the cluster $A$ and the environment $B$ we write the energy functional for the entire system $A+B$ as

[TABLE]

where $\rho(r)=\gamma(r,r)$ is the density, and $E_{H}[\rho]$ the Hartree energy

[TABLE]

with $\mathcal{C}$ denoting the Coulomb operator. The mixed exchange-correlation functional $E^{AB}_{\mathrm{XC}}[\gamma_{A},\gamma_{B}]$ containing both the lower ($B$) and higher ($A$) level functionals can be written as

[TABLE]

To use Eq. (10) in practice within our embedding approach, the contribution due to the interaction between the two subsystems, $E_{\mathrm{XC}}^{\mathrm{int}}$, should be approximated by the lower-level functional ($B$) applied also to the environment, i.e., $E_{\mathrm{XC}}^{\mathrm{int}}[\gamma_{A},\gamma_{B}]=E_{\mathrm{XC}}^{\mathrm{int},B}[\gamma_{A},\gamma_{B}]$ with

[TABLE]

This allows for very expensive functionals to be used in the cluster $A$, including RPA or quantum chemistry approaches. The drawback is that the error introduced in such a mixed approach is difficult to quantify a priori.

For the special case of hybrid functionals chosen in the present work as the “higher-level” functional $A$ the interaction term $E_{\mathrm{XC}}^{\mathrm{int}}$ can be much more accurately approximated by the full hybrid functional itself

[TABLE]

as the most expensive summation $(i,j\in B)$ is constant and hence not relevant for the optimization of orbitals in $A$. The only approximation here is that the orbitals in $B$ are kept frozen. Even when optimizing large supercells, only a subset $N_{A}\times N_{\mathrm{tot}}$ of the full $N_{\mathrm{tot}}\times N_{\mathrm{tot}}$ orbital pairings needs to be calculated to evaluate the relevant exchange contribution $\left<i\;j\right|\mathcal{C}\left|j\;i\right>$, greatly reducing the numerical effort.

The orbitals in $A$ can now be efficiently optimized minimizing the energy functional of Eq. (8), while the orbitals in $B$ are kept frozen. Consequently, during a single optimization, any change in the electronic structure of $A$ due to a more accurate XC functional cannot lead to a redistribution of charge in $B$. The embedding can now be made self-consistent by alternating between subsystems $A$ and $B$ in freeze-and-thaw cycles: after the initial solution of the entire system using the lower-level functional, an orbital rotation is performed to partition into orbital sets $A$ and $B$. Then, starting with $A$, alternatingly one of the subsystems is optimized while the other one is kept frozen. Each cycle that optimizes the orbitals $A$ uses the higher-level XC functional $A$ while each cycle that optimizes the orbitals in $B$ employs the lower-level functional $B$. For the example of polarons in titania discussed below, we find rapid convergence after about six freeze-and-thaw cycles NeugebauerSelfc .

Since the interaction of the cluster region (e.g. defect) with its periodic image needs to be minimized, conventional defect modeling is hampered by the requirement of large supercells. If the bulk material could well be described by conventional XC functionals, and only the defect structure requires more advanced techniques, conventional techniques still require an expensive evaluation of the entire exchange contribution. The embedding procedure outlined above is ideally suited to significantly reduce computational effort while retaining high accuracy.

III Implementation in VASP

We have implemented the embedding scheme outlined above in the Vienna ab initio simulation package (VASP) using the projector augmented wave method of Blöchl in the implementation of Kresse and JoubertBloechl ; Kresse ; KresseFuerthau . Usage is simple: in a first step, a conventional DFT calculation of a system is performed. In a second step, the localized atomic basis functions $\left|\alpha_{i}\right>$ are defined. Our implementation currently supports the PAW basis functions (the pseudo partial waves) and standard spherical harmonics (including hybrid orbitals such as $sp^{3}$) with the radial dependence taken from suitably scaled hydrogen functions. VASP then starts an embedded calculation, performs the orbital rotation and optimizes the set of $N_{\mathrm{emb}}$ orbitals localized on the cluster $A$ while the remaining fully occupied Kohn-Sham orbitals of the environment $B$ are frozen (for spin-polarized calculations the two spin components are handled independently in terms of the rotation and the number of optimized orbitals). During the freeze-and-thaw cycles, no further localization procedure according to Eq. (5) is required, and orbital sets $A$ and $B$ are interchanged. To enforce orthogonality between the currently optimized orbitals and the frozen environment we use the frozen orbitals as projector.

We do not currently support $k$-point sampling in the embedding calculation, since the orbital rotations at different $k$-points are not independent. Likewise, forces are not currently implemented in our formalism. The geometries used in this work were taken from Ref. MerzukDefects, for the defects in silicon and where relaxed using the HSE functional similar to Ref. PolaronCesare, for the polarons in titania.

To compare final energies calculated from the functional Eq. (8), we need to evaluate the full exchange energy with all orbitals once after self-consistently converging the orbitals of $A$. We benchmark our embedding approximation against a fully self consistent optimization of all orbitals using the hybrid functional. Additionally, we compare against evaluating the hybrid functional with orbitals obtained from using a conventional functional (PBE). Obviously, both the embedded and the PBE orbitals are, by construction, not self-consistent with respect to the hybrid functional. In comparison with a full self-consistent optimization, a single evaluation step using the full hybrid functional with non-self consistent orbitals takes a small amount of time, while substantially improving the accuracy: errors in an inexact evaluation of the interaction between subsystems in Eq. (11) are eliminated. We denote corresponding energies by an asterisk (∗) in the following.

IV Results

IV.1 Point defects in silicon

As a first practical test of our new algorithm, we consider point defect structures in silicon MerzukDefects . We use a 64 atom supercell, and neglect $k$-point sampling both in the full hybrid benchmark and in the embedding calculations. Due to the localized nature of the covalent bonds involved, conventional Kohn-Sham DFT fails to correctly reproduce experimental observations. By using hybrid functionals or even more advanced RPA formulations MerzukDefects , these problems are mitigated. However, comparison with more accurate correlation functionals such as RPA and experiment show that currently available methods yield a wide range of predictions depending on the employed functional TkatchenkoExp , highlighting the necessity to move towards higher level correlated wavefunction approaches.

Due to the large supercells required to avoid interaction of the defect sites with their periodic images, embedding the orbitals close to the defect site seems desirable. As benchmark for our embedding scheme, we consider defect formation energies of a set of common interstitial defects and vacancies. We aim to reproduce the energetics of full hybrid functional calculations based on PBEh by a cheaper embedding calculation in which only a few (six to ten) orbitals localized in the immediate vicinity of the defect (taken to be the $A$ orbitals of the cluster) are treated using the hybrid functional, while the remaining  118-122 orbitals (taken to be the $B$ orbitals) are only treated by PBE. We also compare our results to purely DFT-based predictions.

For all Si defect calculations, we choose as atomic orbitals $\left|\alpha_{i}\right>$ the PAW pseudo-partial waves of the Si atoms at and directly adjacent to the defect site, resulting in $N_{A}=16$ (vacancies) or 20 atomic basis functions (1 $s$ + 3 $p$ per atom) for most cases. Increasing the number of basis functions per atom increases the overall overlap of the occupied Kohn-Sham orbitals with the defect site at the cost of a larger $N_{A}$, and thus a larger overlap matrix $W\in\mathbb{C}^{N_{A}\times N}$ of Eq. (3). Consequently, the number of singular values $\sigma_{i}$ increases. To obtain a set of orbitals well localized at the defect site, we choose all orbitals with singular values $\sigma_{i}>0.5$ as embedded orbitals. This procedure yields a number of selected embedded orbitals $N_{\mathrm{emb}}^{A}$ from 9 (X defect) to 16 (T-defect), in line with the number of Si-Si bonds one would expect for the respective defect sites. For example, each of the two defects atoms of the dumpbell defect (X) interacts strongly with four close neighbors in the surrounding lattice and with the other atom in the dumpbell, yielding a total of nine covalent bonds. Indeed, we find nine singular values substantially larger than 0.5 for this defect. The PAW basis functions we choose yield a set of orbitals with a bimodal distribution: a significant number of orbitals with $\sigma_{i}\approx 1$, well seperated from delocalized orbitals with small overlap $\sigma_{i}\approx 0$ with the defect site. The threshold of 0.5 is therefore a good compromise between choosing all $N_{A}$ possible orbitals (which will include orbitals with very small singular values) and too few orbitals that will not allow for reasonable optimization. Note, however, that care must be taken to check that there are no singular values close to 0.5, to avoid arbitrarily including (or discarding) orbitals upon small fluctuations in $\sigma_{i}$.

As mentioned in Sec. II, calculating the defect formation energies from the total energies of Eqs. (6)-(11) for the defects and the defect-free (bulk) system proves challenging. For the embedded case, vacancies and interstitials change the number of orbitals, and make a comparison of absolute energies problematic. We therefore compare the predictions by different methods for the formation energies by evaluating the same hybrid energy functional with the help of all orbitals (i.e., not just the ones localized at the defect) generated by the different methods [denoted by an asterisk (*)].

Our results for various defect structures are summarized in Tab. 1. Overall, we find excellent agreement between the hybrid functional benchmark and our embedding approach. For simple, non-metallic defects such as the dumbbell configuration (X), the hexagonal hollow (H) and a lower-symmetry variant (C3V) we find that both the embedding, as well as the evaluation of the hybrid functional with the low-level DFT orbitals produces good agreement with benchmark calculations (see second column of Tab. 1). By contrast, the metallic tetragonal site is badly described by DFT: it features one interstitial Si atom coordinated to its four nearest neighbors, so that the local coordination of the interstitial is identical to the other Si atoms. This position is unique insofar that the highest occupied orbital is threefold degenerate ($t_{2}$ symmetry) but only occupied by two electrons. This degeneracy is preserved in DFT yielding three fractionally occupied orbitals with occupation numbers $f_{i}=2/3$. Consequently, the evaluation of the hybrid energy functional based on these orbitals fails to yield reasonable formation energies. By contrast, the embedding method locally breaks the degeneracy as does the full hybrid calculation, leading to good agreement of the embedding results with the benchmark (see T, VJT and V lines in Tab. 1).

Our results compare poorly with experimental data: one important reason is the interaction of periodic images of the defects to the supercell size. We therefore consider a larger supercell of 512 atoms, still with a single defect. We note that on our hardware, the full hybrid calculations takes ten times as long as the embedded one, with a relative error of 0.4% in total energy. We find a substantial change in results for the larger cell (compare Tab. 2), that now fit well to experimental results for the H defect. To achieve better agreement also for vacancies (V) requires a more accurate treatment of electronic correlation (e.g., RPA) or inclusion of Van der Waals contributions TkatchenkoExp .

IV.2 Polarons in titania

As a second demonstration of our method, we consider the formation of polarons in titania. We consider a $2\times 2\times 2$ supercell with 24 Ti and 28 O atoms in the rutile structure. In an accurate hybrid functional description, an additional electron localizes, distorting the lattice and forming a small polaron. The distortion decreases the energy compared to a delocalized charge. A full hybrid functional calculation yields a decrease in energy by 514 meV for the distorted geometry. We will use this value in the following as benchmark for our embedded description of the polaron. Concerning the atomic basis functions $\left|\alpha_{i}\right>$ used for the initial localization, we typically use the 1 $s$ and 5 $d$ orbitals of the Ti atom centered at the small polaron deformation, as well as all $s$ and $p$ orbitals of the six nearest neighbor oxygen atoms, yielding a total of $N_{A}=30$ orbitals. Geometries for the distorted structure were relaxed using HSE calculations.

Our results for the small polaron formation energy are summarized in Tab. 3. Conventional density functional theory invoking a PBE functional is not capable of reproducing small polaron formation, predicting even a negative energy gain (i.e. energy costs) of -355 meV to form the polaron. Inserting the DFT orbitals in the hybrid energy functional leads to a correction of the sign. However the energy gain is underestimated by a factor of two (270 meV), see Tab. 3. A single-cycle embedding calculation yields a slightly larger error predicting 190 meV. The origin of this error is obvious: while the hybrid functional tries to localize the charge in the cluster region $A$, the surrounding region $B$ cannot react to the substantial change in the electrostatics, since all $B$ orbitals are frozen.

Subsequent freeze-and-thaw cycles rapidly improve the result: we alternate between optimizing the two sets of orbitals $A$ and $B$, one with the expensive hybrid, the other with pure GGA (PBE). We find convergence in about seven iterations, quite independent of the number of embedded orbitals $N_{\mathrm{emb}}^{A}$ (Tab. 3). As minimum requirement for $N_{\mathrm{emb}}^{A}$, the Ti atom at the center of the distortion and the surrounding oxygen atoms need to be treated accurately, which is already achieved with as few as six orbitals (Tab. 3). Note that only choosing the central Ti atom as atomic basis, $N_{A}=N_{\mathrm{emb}}^{A}=6$, yields a smaller polaron energy than chosing the six orbitals with the highest singular values from the $N_{A}=30$ localized orbitals including also the closest oxygen atoms. The reason is that in the latter case, the response of the surrounding shell of oxygen atoms is - to some degree - also treated by the hybrid functional. However, further increasing the number of localized orbitals $N_{A}$ by, e.g., also including a shell of neighboring Ti atoms does not result in a stronger overlap of $A$ orbitals on the central Ti atom (note that the localization procedure does not distinguish between the different atomic basis functions $\left|\alpha_{i}\right>$). Consequently, such a large $N_{A}$ would require a comparatively large $N_{\mathrm{emb}}^{A}$ to ensure that orbitals close to the central site are included in the embedded calculation. Otherwise accuracy is lost. Indeed, we find a better agreement with the benchmark for $N_{A}=N_{\mathrm{emb}}^{A}=30$ than for $N_{A}=90$, $N_{\mathrm{emb}}^{A}=30$ (see Tab. 3). Since the numerical effort of the embedding calculation scales with a power of $N_{\mathrm{emb}}^{A}$, in practice a small $N_{A}$ that allows for $N_{A}\ll N$ is preferable.

It is intstructive to follow the charge density variations along the freeze-and-thaw cycles [Fig. 1]. Additional charge is localized in the $A$ cycles using the hybrid funcional in the cluster (top row in Fig. 1). The density spreads out again and the environment relaxes in the $B$ cycles when the orbitals $B$ of the environment are optimized using the DFT functional. However, the magnitude of these changes quickly decreases with the iteration number and yields a well-converged density (and well-converged energy) within 7 iterations.

The full hybrid and the converged embedded unpaired spin densities closely match (Fig. 2) (b,c). By contrast, the DFT density does not show a strong localization of the surplus electron at all (Fig. 2) (a). Indeed, projecting the converged polaron orbital (i.e., the occupied majority spin Kohn-Sham orbital with the highest energy) onto the central Ti atom of the distortion yields quite small values for the overlap (0.39) for DFT, while the full hybrid (0.69) and embedded calculations (0.65) agree quite well. This underlines that despite the correct sign for the energy gain when using the DFT orbitals in the hybrid energy functional, the DFT description of the charge density is qualitatively deficient.

V Conclusions

We have demonstrated a new embedding framework based on a suitable rotation in the subspace of fully occupied Kohn-Sham orbitals. Using a projection on local basis functions, a set of orbitals may be localized at a site of interest, for example a defect. Subsequently, these localized orbitals inside the cluster can now be optimized based on a more expensive exchange-correlation functional, such as a hybrid functional involving the exact evaluation of Fock exchange. Since exchange interactions within the frozen environment are neglected, the computation time is drastically reduced. The response of the environment to the charge rearrangement in the cluster can be self-consistently included by freeze-thaw cycles in which alternatingly the orbitals in the embedded cluster or in the environment are optimized.

We have implemented our ansatz in the popular VASP software package. As proof of principle, we have applied our method to two problems of interest: a set of defects in bulk silicon, and small polarons in bulk titania. We find excellent agreement with (much more expensive) benchmark bulk hybrid calculations.

Acknowledgements.

The authors gratefully acknowledge support by the FWF via the SFB-41 ViCoM.