Implementation of density functional embedding theory within the projector-augmented-wave method and applications to semiconductor defect states

Point defects universally exist in all types of solid state materials and play important roles in engineering material properties. The presence of defects significantly affects the chemical, electrical, optical, and mechanical properties of the material.¹ Understanding the physics of point defects is thus crucial for designing new materials for various applications. In particular, defects have been shown to be one of the key factors determining the performance of thin film solar cells:² on the one hand, certain levels of intrinsic defects or extrinsic dopants are the source of free carriers, increasing the conductivity of the semiconductor; on the other hand, wrong types of defects also serve as non-radiative recombination centers and introduce immobile trap states. It is thus of paramount importance to understand defect states in order to ultimately control them. The dearth of experimental techniques to probe bulk atomic-scale defects motivates use of first-principles quantum mechanics calculations to do so.

Density functional theory (DFT) based on plane-wave (PW) basis sets is currently the state-of-the-art in the material science community for describing periodic extended systems. Although being successful in general, DFT suffers from approximate treatment of electron exchange-correlation (XC), resulting in self-interaction error³ and the related underestimation of band gaps. In quantum chemistry, more accurate correlated wavefunction (CW) methods (Møller-Plesset (MP) perturbation theory, configuration interaction (CI), coupled cluster (CC), etc.) are available, but these methods are not straightforward to implement in a periodic setting and much more expensive within the PW formalism.⁴ This is especially true for the simulation of isolated defects, which usually requires large supercells to avoid interactions between periodic images. The poor scaling behavior of the PW-CW methods significantly limits their applicability in these scenarios. This problem becomes even more severe for charged defects, which requires careful consideration using periodic boundary conditions (PBCs). For example, in some solar cell materials, the charges carried by the defects create local electrostatic potential fluctuations, which can be a major reason for low open-circuit voltages.⁵ To study such charged defects, besides large supercells, specially designed corrections^6–9 are also needed to diminish the long-range electrostatic interactions between periodic images.

A promising alternative to eschew these difficulties associated with the periodic PW formalism is the use of embedded cluster calculations. In this approach, we exploit the locality of the defect and construct a finite cluster model, representing the critical region of interest. The influence of the surrounding parts of the system (the environment) is represented by an embedding potential. The environment and its interaction with the cluster can be computed using a lower level method such as a DFT approximation, while the cluster can be studied using high-level embedded correlated wavefunction (ECW) methods. Various approaches to constructing the embedding potential have been proposed, which have been summarized elsewhere.^10–12 In this work, we focus on the density functional embedding scheme developed by Huang et al.¹³ In this scheme, one can define a unique embedding potential for both subsystems, which exactly reproduces the electron density of the total system by summation of the electron densities of the two embedded subsystems. By solving an optimized effective potential (OEP) problem, we obtain this unique embedding potential at the DFT level (see below for details). We can now conduct ECW calculations on the embedded cluster without PBCs, using the well-developed formalisms of quantum chemistry. The density functional embedding scheme has been successfully applied to studying absorbates (such as H₂ or O₂) on metal surfaces, including Al and Au.^14–16 While it has been shown that the embedded metal cluster can decently reproduce bulk metal properties, the accuracy of the density functional embedding method beyond metallic systems has been tested sparingly.¹⁷ Several other embedding schemes have been applied to systems featuring more covalent bonding character, in both molecules and condensed matter.^{10–12,18–22} We therefore would also like to extend the scope of the density functional embedding theory (DFET) method beyond metals. In this paper, we will demonstrate some preliminary test cases illustrating the potential capability of the density functional embedding method for the study of semiconductor defect systems.

While performing PW-DFT calculations on semiconductors, there are two primary ways to handle the large wavefunction amplitudes near the core region that otherwise would require very large kinetic energy cutoffs using PWs: pseudopotentials (PPs) and the projector-augmented-wave (PAW)^23,24 method. The density functional embedding method was implemented previously in the ABINIT^25–27 program, in conjunction with norm-conserving PPs.²⁸ Meanwhile, another highly parallelized program that is widely utilized in the material science community is VASP.^29–32 VASP provides a wide selection of XC functionals, ranging from local density approximation (LDA)³³ to hybrid functionals such as Heyd-Scuseria-Ernzerhof 06 (HSE06).³⁴ Compared to ABINIT, VASP also provides a simpler U-J correction formalism developed by Dudarev et al.,³⁵ which is often employed to alleviate the self-interaction errors for transition metal elements in semiconductor and ionic materials. More importantly, VASP offers a robust, well-tested library of PAW potentials that covers the entire periodic table. The PAW formalism is more accurate than PPs, as it (i) offers an exact transformation between the original, oscillating wavefunction and the pseudized, slowly varying pseudo-wavefunction used in the computation and (ii) explicitly treats the all-electron wavefunctions, albeit within the frozen-core approximation. The implementation of PAW in VASP offers an efficient way to conduct frozen-core, all-electron calculations with reasonably low kinetic energy cutoffs (typically even lower than the norm-conserving PP calculations). Considering these advantages, we would like to extend our implementation of the density functional embedding scheme to VASP. This extension should enable us to conduct embedding calculations with the PAW-DFT formalism beyond the current PP setup. As we will show, the PAW formalism introduces additional complications into the embedding scheme, which require a series of modifications of the embedding algorithm. In the following, we describe the algorithm developed in the new implementation and we test the robustness of the new algorithm on various types of systems, including covalently bound molecules, metals, and most importantly, semiconductors including defects.

In the following, we briefly review the algorithm and the numerical details of the density functional embedding method here. Interested readers can find a detailed treatment in Ref. 13 for more complete theoretical derivations.

Consider a total system partitioned into two subsystems, labeled as system A and system B, respectively. The essential idea of the DFET is to find a unique embedding potential V_emb for both systems A and B, such that

The $n K V emb ,K=A,B$⁠, represent the ground-state electron densities of the two subsystems obtained from separate, converged self-consistent field (SCF) calculations in the presence of V_emb and n_ref represents the ground-state electron density of the total system, computed without embedding potential. In this work, we always use the same embedding potential for different spin states, and consistently we always perform non-spin-polarized calculations to maintain the spin symmetry. In practice, in order to obtain V_emb, we need to maximize the following extended Wu-Yang (WY) functional³⁶ with respect to V_emb:

Here, E_K[n_K], K = A, B is the energy functional of system K, including its interaction with $V emb ( ∫ V emb ( r → ) n K d r → )$⁠. Given that

it can be proven that the derivative of the WY functional with respect to V_emb is simply¹³ 

Equation (4) gives the gradient used in the optimization procedure: when W reaches its maximum, the gradient vanishes and condition (1) is naturally satisfied. The entire procedure is effectively solving for an OEP, which, in our case, turns to be the optimal unique embedding potential.

Analytically, Eq. (3) is guaranteed by the variational principle. In the current implementation in ABINIT, all the functions (densities, potentials, etc.) are expressed on one uniform three-dimensional grid and all integrations are performed as summations over this grid. In this circumstance, Eq. (3) always holds numerically, even using a finite grid size. However, as we discuss below, this is no longer true for dual grid schemes.

One of the key steps of the density functional embedding procedure is to perform SCF-DFT calculations for both subsystems embedded in an external potential V_emb, from which both densities (n_K) and energies (E_K[n_K]) are obtained. Therefore, we need to consider how to conduct PAW-DFT calculations in conjunction with the presence of an arbitrary external potential. Arbitrary in this context refers to the typical length scale of potential fluctuations: while including a slowly varying potential (on the length scale of interatomic distances) is quite straightforward, rapid potential fluctuations close to the ionic cores affect the PAW projection and thus need to be treated carefully. In this section, we briefly review the PAW formalism and discuss the necessary modifications introduced by the arbitrary external potential. The readers are referred to previous papers^23,24 for a more detailed description of the PAW method.

Using Blöchl’s nomenclature,¹⁸ the exact all-electron (AE) wavefunction of the n-th band $Ψ n$ is expressed as a soft nodeless pseudo (PS)-wavefunction $Ψ ̃ n$ via Eq. (5),

where $ϕ i$ and $ϕ ̃ i$ are the onsite AE and PS basis sets, respectively. These two sets of basis functions are identical beyond the core radius R_aug but differ inside the core region. $p ̃ i$ are the projector functions that are orthonormal to the PS basis: $p ̃ i ϕ ̃ j = δ i j$⁠. The projection introduced by Eq. (5) damps the sharp oscillations of the wavefunction within the core radius R_aug such that the resulting soft PS wavefunction $Ψ ̃ n$ can be expanded using considerably fewer plane waves. Usually, the frozen-core approximation is employed, i.e., the AE wavefunctions of the core electrons are fixed during the SCF iterations. In contrast to norm-conserving PPs, the exact AE wavefunctions for both valence and core electrons are explicitly represented and expanded using the ansatz given by Eq. (5). Therefore, the PAW scheme can effectively be considered as a frozen-core, all-electron level of theory.

In correspondence with the decomposition of the wavefunction, the total density of the system is decomposed into three parts,

The first term $n ̃$ is the PS density calculated directly from PS wavefunctions on the uniform PW grid. The second (n¹) and third (⁠$n ̃ 1$⁠) terms are the one-center densities (denoted by the superscript 1) defined within the augmentation sphere, associated with the AE and PS basis functions, respectively. Due to the rapidly oscillating character of the one-center terms, both densities are computed on dense radial grids centered on the ion positions.

Within the PAW ansatz, the Kohn-Sham (KS) Hamiltonian of the system can be written as

Here, $v ̃ eff$ is the Hartree and XC potential due to the soft PS valence electron density $n ̃$⁠, the PS core electron density plus the nuclear charges $n ̃ Z c = n ̃ c + n Z$⁠, and the compensation charge $n ˆ$⁠. All these terms are soft (meaning they have no sharp oscillations) so they can be evaluated and stored on a coarse, uniform grid. The compensation charge $n ˆ$ and its angularly decomposed components $Q ˆ i j L ( r → )$ are introduced to compensate the different multipole moments of the AE and PS densities within the core region.²⁴ The term $v ̃ eff 1$ in Eq. (10) is the one-center counterpart of $v ̃ eff$⁠, which is still soft but evaluated on much finer, ion-centered radial grids. The sharp one-center potential $v eff 1$ corresponds to the oscillating components of both valence and nuclear-core charges (n¹ and n_Zc). Naturally, these strongly varying terms are computed and stored on the fine, ion-centered radial grids.

When we introduce an embedding potential in the system, it enters all three potential expressions, resulting in extra terms (compare to Eqs. (43), (45), and (46) of Ref. 24),

where $n ̃ c$ and n_c are the PS and AE core electron densities and $n ̃ Z c$ and n_Zc are the core electron densities plus the nuclear charges. It is straightforward to perform the addition in Eq. (11), as both the original soft potential $v ̃ eff$ and the embedding potential V_emb are stored on the same uniform grid. However, in Eqs. (12) and (13), the one-center potentials (⁠$v eff 1$ and $v ̃ eff 1$⁠) are expressed on an angularly decomposed radial grid. Therefore, in order to incorporate the embedding potential, an extra step is needed to transform V_emb from the uniform grid to radial grid, using the following projection:

In this equation, S_LM are the real spherical harmonics, $V emb LM$ is the LM-component of the embedding potential, in which L and M denote the angular momentum of the corresponding channel. r and $r ˆ$ are the length and direction of the vector $r → − R → I$⁠, and $R → I$ is the ion position. Then, the angularly decomposed $V emb LM$ is added in Eqs. (12) and (13) to account for the effects of the embedding potential on the one-center terms in PAW.

The dual grid transformation in Eq. (14) has profound consequences for the evaluation of the energy derivatives with respect to the embedding potential. According to Ref. 24, the total KS energy is

Here, the first term is the summation of the orbital eigenvalues, which double counts the electron-electron interactions. In order to correct for this error, we have to include the double counting terms, which are denoted with subscript “dc.” The total energy also includes the ion-ion interaction term, which is represented by $U ( R ⃗ I )$⁠. Considering the variational principle and using the definition of $n ˆ$ (Eq. (27) in Ref. 24), the derivative of the total energy is

in which the onsite density matrix ρ_ij is defined as

Here, the double counting terms (denoted with subscript “dc”) and the ion-ion interaction term $U ( R ⃗ I )$ vanish as they do not contain V_emb explicitly. We only consider the valence bands for the summation over n since we only aim to match the valence electron densities, as usual within the frozen-core approximation. The first two terms of Eq. (16) are the soft charges given by the PS wavefunction, which are easy to compute. However, the exact evaluation of the third term is less straightforward, as V_emb enters $v eff 1$ and $v ̃ eff 1$ in an indirect way. In consistent with the radial grid used in the VASP code, we insert Eqs. (9), (10), (12), and (13) into Eq. (16) and evaluate all the integrals using radial coordinates. In this way, we reach the following expression for the third term in Eq. (16):

where i = (l, m) and j = (l′, m′) denote the angular momenta of the projectors (⁠$p ̃ i$ and $p ̃ j$⁠). $C l m l ′ m ′ LM$ is the integration of the product of three real spherical harmonics (⁠$C l m l ′ m ′ LM =∫ S LM ( r ˆ ) S l m ( r ˆ ) S l ′ m ′ ( r ˆ ) d r ˆ$⁠), which can be easily calculated via linear combinations of Clebsch-Gordan coefficients. L,  M are the angular indices for the angular components of the embedding potential and Q_lml′m′ and $Q l m l ′ m ′ L$ are the depletion charges,

The definition of the compensation function g_L(r) is given by Eq. (61) in Ref. 24. In Eq. (18), the most critical step is evaluating $∂ V emb LM ( r ) ∂ V emb ( r → )$⁠, which defines how the values on the radial grid vary with respect to the values on the uniform grid. The speed of computing this term depends on the projection algorithm we use for Eq. (14). We will compare different projection algorithms in more detail in Secs. II C and II D after we clarify the modifications we need to make to the OEP procedure.

As can be proven analytically, the energy derivative given by Eq. (16) should be identical to the AE density $n ( r → )$⁠, consistent with Eq. (3). However, at a typical PW grid density, the energy derivative significantly differs from the AE density given on the uniform grid. In Fig. 1, we demonstrate this difference using the Cl₂ molecule as an example (for more computational details, see Sec. III). It is clear that the AE density $n ( r → )$ exactly corresponds to the energy derivative in the interstitial region but shows much stronger variations in the core region.

The reason for this discrepancy is directly related to the dual grid scheme adopted in the PAW algorithm. Since the radial grid is much finer than the uniform grid, the transformation from the uniform grid to the radial grid must use an interpolation. Consequently, the change of a potential value at one uniform grid point effectively leads to changes of potential values at multiple radial grid points. In other words, even though the derivative $∂ V emb LM ( r ) ∂ V emb ( r → )$ should have rigorously a δ-function form (⁠$F ( r ) δ r − r →$⁠), it is slightly nonlocal in practice. Therefore, the resulting energy change is no longer determined by the density value on one particular point (⁠$n ( r → )$⁠), but rather a smeared average of the densities of the neighboring region, explaining why the energy derivative grid is always softer than the AE density grid.

As a consequence, Eq. (3) breaks down and the WY functional given by Eq. (2) is no longer appropriate, as the simple formula for the derivative given by Eq. (4) no longer holds. In practice, the difference between the densities and the derivatives is so large that no meaningful optimizations can be done if we use densities as an approximate derivative. To resolve this problem, we slightly modify the WY functional,

With the associated derivative,

Basically, we replace the densities with the energy derivatives. Instead of matching the densities, we match the exact energy derivatives, which is physically similar since the exact energy derivative is effectively a smeared average of the density. We will show in Sec. III that this modification indeed does not change the essential physics of the original density functional embedding method. To conduct the optimization, we have to compute the exact energy derivatives for both the total system and the subsystems in each iteration, using Eqs. (16) and (18).

As we discussed in Sec. II B, we need to transform the embedding potential from the uniform PW grid to the PAW radial grids (Eq. (14)). This projection can be performed using a reciprocal-space algorithm, which starts with the Fourier interpolation of the embedding potential,

Then, we utilize³⁷ 

to expand the phase factor $e i q → k ⋅ ( r → − R → I )$ in spherical harmonics and spherical Bessel functions, obtaining

Here, we define $r ˆ$ as the direction of $r → − R → I$⁠. Comparing it with Eq. (14) yields

Here, $r= r → − R → I$⁠, and $r → n$ represents the n-th uniform grid point in real space, with n = (i, j, k) ranging from (1,1,1) to (N_x, N_y, N_z). $R → I$ is the coordinate of the ion position. N is the total number of uniform grid points (N = N_xN_yN_z), and $q → k$ is the k-th point in reciprocal space. With this projection formalism, we can easily derive an expression for

The reciprocal-space projection algorithm is a very robust algorithm that we found well-behaved in all scenarios. However, the derivative given by Eq. (26) is extraordinarily expensive to compute. The summation in Eq. (26) must be conducted for every radial grid point r and every uniform grid point $r → n$ near every ion $R → I$⁠. The summation over the entire reciprocal space (⁠$∑ q → k$⁠) contains millions of terms, making the treatment of realistic application cases challenging. Nevertheless, the robustness of the reciprocal-space projection algorithm makes it a perfect benchmark method for more approximate projection algorithms. Therefore, we used this algorithm to compute the exact derivatives in Fig. 1.

The intrinsic reason for the heavy computational cost of the reciprocal-space projection algorithm is the nonlocality of the Fourier interpolation given by Eq. (22). Basically, the potential value $V ( r → )$ at an arbitrary position depends on the potential values on all the reciprocal-space grid points (⁠$v ̃ q → k$⁠), which are again determined by the values on all the real-space grid points (⁠$V ( r → n )$⁠). In other words, we need the information of the entire potential grid before we can find the interpolated potential value at a particular position. This is counter-intuitive, as $V ( r → )$ should only be affected by the $V ( r → n )$ values on the grid points $r → n$ that are spatially close to $r →$⁠. To overcome the nonlocality associated with the Fourier transformation, we need a projection algorithm based on a real-space interpolation scheme.

The algorithm we adopt is a modified version of Misner’s decomposition method.³⁸ The basic idea is to first expand the potential using both angular and radial basis sets,

Then, we reassemble the radial expansion to obtain $V emb LM r$⁠,

Here, R_K(r) are orthonormal radial basis functions (vide infra) and V_LMK are expansion coefficients computed via integrations on a uniform grid,

To conduct the integration in Eq. (29), we introduce a secondary uniform grid ${ r → m }$ for each ionic core, which is only defined within the augmentation sphere Ω defined by R_aug. The introduction of the secondary grid is the first major improvement we made to Misner’s original method. It is necessary as the primary PW grid ${ r → n }$ is usually not fine enough for the evaluation of V_LMK. In this work, we always adopt a comparatively dense dimension of 50 × 50 × 50 for the secondary grid. To interpolate the potential values on the secondary grid (⁠$V ( r → m )$⁠) from the primary grid (⁠$V ( r → n )$⁠), we utilize the uniform cubic B-spline scheme³⁹ generalized to 3D space. Basically, we use the localized 1D cubic B-spline functions defined for each primary grid point as the fundamental basis set,

We assume that our 3D continuous potential can be expanded using these basis functions,

Once we have the expansion coefficients f_ijk, we can compute the potential values on the secondary grid $r → m$ using Eq. (31). The coefficients f_ijk are chosen such that the continuous interpolated $V ( r → )$ matches all the values on the primary grid $V ( r → n )$⁠. In practice, we need to construct a series of 1D δ-grids for each of the dimensions to obtain f_ijk,

Following the standard 1D interpolation procedure,³⁹ we interpolate the i-th 1D δ-grid in the x direction by solving the following equation for the 1D interpolation coefficients ${ a i ′ i }$⁠:

Note that due to permutation symmetry, Eq. (33) needs to be solved only once for all the i and i′. The resulting coefficients are very localized, given that $a i ′ i$ vanishes quickly with respect to increasing $i − i ′$⁠. Similarly, we can obtain the interpolation coefficients ${ b j ′ j }$ and ${ c k ′ k }$ for y and z directions, respectively. Then, we construct the interpolation coefficients for the 3D δ-grids,

The final coefficients f_ijk for $V ( r → )$ are computed using Eq. (35), recognizing that the $V ( r → )$ grid can be written as a linear combination of the 3D δ-grids, $V ( r → ) = ∑ i j k V ( r → i j k ) δ i j k ( r → )$⁠,

With this procedure, it is also straightforward to evaluate the term $∂ V emb LM ( r ) ∂ V emb ( r → )$ by exploiting the linearity of the interpolation scheme. When we consider $∂ V emb LM ( r ) ∂ V emb ( r → n )$ for a particular uniform grid point n = (i, j, k), we only need to replace the term $V ( r → m )$ in Eq. (29) with the δ-grid $δ i j k ( r → m )$ defined in Eq. (34) and insert it into Eq. (28). The resulting radial function from Eq. (28) is exactly

Formally, it is not obvious why this procedure would be more efficient in evaluating derivatives, as the second summation in Eq. (36) still involves 50 × 50 × 50 terms, i.e., the size of the secondary grid. However, noticing that the interpolated function for the 3D δ-grids (⁠$δ i j k ( r → )$⁠) is strictly localized around the primary grid point $r → n = ( ijk )$⁠, the summation in Eq. (36) can be substantially reduced. With simple distance cutoffs, most of the terms with large $r → m − r → n$ can be dropped without losing accuracy. In fact, the locality of the interpolated δ-grid is a reflection of the locality of the real-space cubic B-spline interpolation algorithm, which is the intrinsic reason for the orders of magnitude speedup provided by the real-space projection procedure.

Note that $∂ V emb LM ( r ) ∂ V emb ( r → n )$⁠, as given by Eq. (36), and thus Eq. (18), is independent of $V ( r → )$⁠, further reducing the computational cost of the interpolation scheme. As long as the grid setup, the PAW library setup, and the atomic geometry are unchanged, the results from Eq. (18) remain constant. Consequently, during the OEP optimization, Eq. (18) needs to be computed only once in the beginning, and the results can be recycled in further iterations. This property further speeds up the code by several orders of magnitude.

In this section, we introduce several important aspects about the implementation of the real-space projection algorithm.

The most important factor in the real-space projection is the selection of the radial basis sets R_K(r). Instead of the r-scaled Legendre polynomials used in the original Misner method, we adopt sinc functions,

where sinc functions are orthonormal within [0, R_cut], which is a fundamental requirement for R_K(r),

In the r → 0 limit, all the sinc functions have well-defined finite limits, in contrast to the r-scaled Legendre polynomials, which diverge at short range. Because of this property of the sinc functions, the integration in Eq. (29) is well-defined within the entire augmentation sphere. By contrast, in the original Misner method, the real-space integration can only be performed within a shell to avoid the numerical instability at r → 0. This is an another major improvement we made to the original Misner method.

Although the above procedure is theoretically feasible, we find that the resulting $V emb LM$ is still numerically inaccurate at very short range. Increasing the number of radial basis functions (K_max) improves the accuracy but slows down the computation. Furthermore, higher K values introduce strongly oscillating basis functions, which require a much finer secondary grid for the real-space integration. Compared to a radial grid, the uniform Cartesian grid lacks the necessary accuracy close to the center of the sphere, introducing larger numerical errors. This is an intrinsic problem for uniform grid layouts and should be avoided also from other considerations. If we Taylor-expand the potential and project it onto different angular channels, we can derive the short-range behavior of $V emb LM$ analytically,

To avoid the problem stated above, we switch to Eq. (39) in the limit r < ϵ. The value ϵ is usually inversely proportional to K_max. As a rule of thumb, we take ϵ = R_cut/K_max, with K_max set to 20 in this paper, which we found high enough for our applications. The prefactor α is chosen to enforce continuity of $V emb LM ( r )$⁠.

In the limit r → R_cut, all sinc functions rigorously converge to zero and thus can only expand functions with the same boundary condition. In general, our arbitrary embedding potential has finite values over the entire space and thus cannot be expanded by the sinc functions in a naïve way. To solve this problem, we introduce a padding region outside of the augmentation spheres, within which the potential is gradually switched off by a cosine function,

Here, d is the thickness of the padding region, which equals to 0.1 Å in this work. Instead of $V emb ( r → )$⁠, we actually project $V ̃ ( r → )$ and define the cutoff radius as R_cut = R_aug + d. After the projection, we discard the part of $V emb LM ( r )$ in the padding region and only keep the part within R_aug in the PAW equations.

In Eq. (16), the soft charges (⁠$n ̃ + n ˆ$⁠) are almost negligible in computational cost, while the third term which corresponds to the core correction is much more expensive. Although the third term exists only in the augmentation sphere, it is still challenging to compute the corrections for all the points inside the sphere. Usually, at the boundary of the augmentation sphere, the AE and PS wavefunctions are already very close, so the derivative correction is small. Therefore, we introduce a scaling factor F_corr. The derivative corrections are only computed for points within F_corr ⋅ R_aug. In practice, F_corr controls the balance between accuracy and computational cost: for larger F_corr, the derivative is more accurate and the calculation is more expensive. We have found a reasonable selection of F_corr to usually lie in 0.67–0.90, with no significant effects on final results within this range.

Another problem we need to resolve during the OEP optimization is the unphysical oscillation of the resulting embedding potential. This problem is well-known for OEP optimizations⁴⁰ and caused by the unbalanced basis sets for the density and the potential. Furthermore, the frozen core orbitals cannot react to the potential within the PAW sphere, possibly leading to unphysically large potential fluctuations. Using the formula in Ref. 40, we introduce a penalty function in the WY functional, damping potential oscillations,

Mathematically, the penalty function unfortunately damps both physical and unphysical potential fluctuations uniformly. Therefore, the contribution of the penalty function (controlled by the parameter λ) should not be too large, lest physical potential oscillations are damped out. In this work, unless stated otherwise, we always set λ to 10⁻⁵ eV/(V²/Ų), which we have empirically found to give robust results. We will show that the physical results derived from the embedding potential are relatively insensitive to the exact value of λ. Nevertheless, a more rigorous approach will be needed to solve the ill-posed OEP problem in PW basis sets, similar to algorithms put forward for Gaussian basis sets.⁴¹ 

In this section, we discuss the algorithms we implemented for the embedded cluster calculations using all-electron GTO basis sets in NWCHEM.⁴² Considering a system embedded in an arbitrary external potential V, we only have to modify the one-electron integrals to include the effects of V in a quantum chemical calculation. The new terms introduced by V are

Here, G₁ and G₂ are the two primitive Gaussians centered at atom positions $R → I 1$ and $R → I 2$⁠, with exponents of a₁ and a₂, respectively.

Normally, this integral can be directly evaluated via a summation over the uniform grid,

However, this real-space summation is problematic if the term $G 1 ( r → ) G 2 ( r → )$ is sharp in the corresponding region. Particularly, when performing all-electron calculations, the basis functions describing the core electrons are usually associated with very large exponents. The typical grid densities adopted in the VASP-PAW setup are certainly not fine enough for computing Eq. (43). Notice that $I∝∫exp − ( a 1 + a 2 ) ( r → − a 1 R → I 1 + a 2 R → I 2 a 1 + a 2 ) V ( r → ) d r →$⁠, so the sharpness of the integrand is directly controlled by the total exponent a₁ + a₂. Therefore, we adopt an alternative reciprocal-space scheme if a₁ + a₂ is larger than 20 a.u. and $R → I 1 = R → I 2 = R → I$ (note that when the Gaussians are sharp and $R → I 1 ≠ R → I 2$⁠, they do not overlap and the integrand vanishes, so we only consider the $R → I 1 = R → I 2$ situation here).

We assume that Cartesian functions are used for the angular part of the primitive Gaussians. So, the two primitive Gaussians can be written as

We perform a Fourier transformation of $V ( r → )$ (note that the potential is relatively smooth compared to the density, so no ultra-fine uniform grid is needed for this Fourier transform) and plug it into Eq. (42). Then, we reach the following final working equation:

Here, $v ̃ ( k → )$ is the Fourier transform of $V ( r → )$ and $k →$ represents all the reciprocal-space vectors. The integrals $f k , l , a =∫ e i k x x l e − a x 2 dx$ can be derived using the following recursion relations:

In this work, we always use Cartesian functions for the angular part of the basis sets. Theoretically, when spherical functions are used, we can still perform the calculations in Cartesians first and transform them to spherical form.⁴³ Nevertheless, once we obtain the modified one-electron integrals, we conduct the rest of the calculations using the existing code. We are thus able to carry out embedded quantum chemistry calculations using any software with the modified one-electron integrals.

To sum up, with the real-space projection algorithm, we are able to evaluate the exact energy derivatives (⁠$δ E δ V ( r → )$⁠) of an embedded system in a very efficient way. With the energy derivatives of the two subsystems and the total reference system at hand, we can maximize the WY functional given by Eq. (20). This maximization eventually leads to the density functional embedding potential at the PAW level of theory. All algorithms mentioned in this section were implemented in a modified version of VASP 5.3.3.^29–32 For the embedded cluster calculations using GTO basis sets, we also developed an accurate reciprocal-space algorithm to compute the effect of an arbitrary embedding potential within the one-electron integrals. This algorithm has been implemented in NWCHEM-6.5,⁴² enabling us to conduct embedded cluster calculations using all-electron GTO basis sets. In Sec. III, we use the code to compute the density functional embedding potentials for various test systems; we also demonstrate the robustness of the algorithm and the embedding method.

In this section, we present the results of several test cases. Unless stated otherwise, all calculations were conducted using a modified version of VASP with the Perdew-Burke-Ernzerhof (PBE)⁴⁴ XC functional and the PAW method. All calculations were conducted using only Γ-point sampling in reciprocal space as we only study molecules and large supercells. The readers are referred to the supplementary material⁴⁵ for more computational details for each test case.

The first test case we present is the covalently bound Cl₂ molecule, serving as proof-of-principle for our approach. Two chlorine atoms are placed parallel to the z-axis, at the center of a 10 Å × 10 Å × 10 Å box, with a bond length of 2 Å. Naturally, the two subsystems are the two chlorine atoms, respectively.

We first examine the quality of the real-space projection using the reciprocal-space projection results as a benchmark. We take the optimized embedding potential of Cl₂ and decompose it into different angular momentum channels using either real-space or reciprocal-space projection algorithms (see Fig. 2).

At short range, the potential is sharper and thus more sensitive to details of the interpolation, so small deviations can be observed between the two algorithms. Nevertheless, overall, the real-space algorithm reproduces the more accurate reciprocal-space results with very high fidelity and much lower computational cost. This confirms the numerical accuracy of the algorithm and validates our selection of projection parameters introduced in Sec. II.

To investigate the effect of the penalty function, we compare converged embedding potentials with and without applying a penalty function (see Fig. 3). Without a penalty function, the resulting potential features unphysical oscillations. The penalty function smoothens the potential, yielding a much more physical shape, with peaks (electron depletion) located in the core regions and valleys (electron concentrations) in the σ-bond region.

More systematically, we conducted the OEP optimization with different λ values and investigated the change of the optimized WY functional values with respect to λ (see Fig. S1⁴⁵). We observed a sharp decrease of the value of the WY functional below λ = 10⁻⁵ eV/(V²/Ų), mainly corresponding to the damping of unphysical oscillations. As we discussed above (Fig. 3), at λ ≈ 10⁻⁵ eV/(V²/Ų), a smooth potential can be obtained, while larger λ values start to damp the physical variations, which features a slower increase of the WY functional value. Therefore, we will adopt λ = 10⁻⁵ eV/(V²/Ų) in the following study.

As discussed above, care must be taken not to remove physical oscillations: the figure of merit here is the final correspondence between the sum of the subsystem densities and the reference density and the transferability of the final embedding potential to other approaches (i.e., quantum chemistry codes). For the OEP application required for density functional embedding, we can thus easily assess the quality of the embedding potential via comparing the densities/energy derivatives of the subsystems to the reference system (see Fig. 4).

Before optimization (i.e., with V_emb = 0), significant differences between the subsystem densities and the reference system densities appear (the root-mean-square deviation (RMSD) equals 1.1 × 10⁻² e/ų for the energy derivative and 1.0 × 10⁻² e/ų for the soft charge). This is expected as the interactions between the two subsystems are completely missing without embedding potential. After the optimization, the energy derivative errors are significantly reduced (RMSD = 3.4 × 10⁻⁴ e/ų without penalty function), which is expected as we were explicitly optimizing them. Moreover, the differences in densities on the real-space grids are also greatly reduced (RMSD = 4.6 × 10⁻⁴ e/ų without penalty function), just as required by the original definition of the density functional embedding scheme. This illustrates that our modifications to the original scheme defined by Eq. (20) do not change the essential physical idea of the original method.

We also note that the penalty function introduces small but non-negligible residual errors in the core region for δE/δV. For the entire 3D grid, the RMSD of the energy derivatives increases from 3.4 × 10⁻⁴ e/ų to 7.1 × 10⁻⁴ e/ų when using a penalty function. By introducing the penalty function, we effectively introduce some extra constraints on the potential, and thus, the density fitting is numerically less accurate. This shows that the penalty function not only damps unphysical but also a certain level of physical variations of the embedding potential, which is almost inevitable. However, the penalty function is necessary to remove the unphysical oscillations, which is critical for the robustness and the transferability of the resulting potential. In order to obtain a potential that can be utilized in conjunction with different basis sets and methods, we have to filter out the oscillations on the length scale of the uniform grid size. As we showed in Fig. 3, the penalty function accomplished this goal, though with a little compromise in reproducing the reference densities.

We conduct the density functional embedding calculations for an Al₁₂ cluster in a 5 × 5 × 4 Al slab surface (see Figs. 4(a) and 4(b)), which was previously studied using our ABINIT implementation.¹⁴ Using this well-understood metal system, we will verify that our new implementation is able to reproduce the same physics demonstrated by our previous implementations. We would also like to understand the effects of the PAW scheme on the numerical details of the resulting embedding potential, compared to the PP calculations. To be consistent with our previous work,¹⁴ no penalty functions are applied in this case. The embedding potentials obtained from both VASP-PAW and ABINIT-PP implementations are shown in Figs. 5(c)-5(e).

The ABINIT-PP potential (Fig. 5(c)) features negative values mainly in the boundary region between the cluster and the environment, representing the attractive interaction due to the metal bonds. Comparison with the VASP-PAW potential (Fig. 5(d)) reveals two distinct regions: in the core region, the potential obtained from the PAW calculations shows many more features compared to the potential from the PP calculations, as expected since the all-electron PAW approach yields a more complete description of the core region. In the interstitial region, the VASP-PAW potential and the ABINIT-PP potential show excellent agreement, as their isosurfaces are similar in both size and shape. For a more systematic and quantitative comparison for the interstitial region, we make a scatter plot comparing the two embedding potentials on a point-by-point basis (Fig. 5(e)). Due to the effects of PAW scheme on the interstitial regions, we find a root mean square deviation of 0.2 V between the two potentials. The PAW wavefunction gives slightly different densities and response properties compared to the PP wavefunction, so some numerical differences are expected. Furthermore, we can never completely rule out the effects of unphysical oscillations, which possibly vary with respect to different core region treatment methods. Nevertheless, the nice overall correlation (R = 0.96) between the two embedding potentials demonstrates that both potentials capture the same physics in the interstitial region, consistent with the isosurface plots.

In summary, comparing embedding potentials from the VASP-PAW and ABINIT-PP calculations shows that our new implementation is consistent with the old PP-based implementation. Moreover, we are able to obtain new features in the core region representing the embedding effects on the core electrons, which were missing from the previous PP calculations.

In this section, we examine the capabilities of the density functional embedding method for treating defect states in semiconductors. We use a 2 × 2 × 2 64-atom supercell of zincblende ZnS as our test system and introduce one Sn_Zn (Cu_Zn) antisite in the center to represent an n-type (p-type) point defect (see Fig. 6). In this case, we always apply a penalty function to the WY functional, so we can obtain more physical and more robust potentials.

For the Sn-doped system, we carve out a cluster containing the central Sn defect and its first and second nearest neighbors, i.e., a SnZn₁₂S₄ cluster. For the Cu-doped system, we adopt the same cluster size, generating a CuZn₁₂S₄ cluster. In the previous cases (homogeneous diatomic molecules and metals), it was natural to assume that both subsystems were neutral. However, the partitioning of electrons is much less straightforward in semiconductors: if we assume covalent bonds between Zn and S, it is reasonable to cleave the Zn–S bond homogeneously; if we emphasize the ionic character of the Zn–S bonds, we should assign all the bonding electrons to sulfur. It is difficult to determine which way is more physical as the Zn–S bonds have both covalent and ionic characters. The DFET requires an appropriate partition as an input and, theoretically, the optimization process should work for any electron or atom partitions. Therefore, in the framework of the density functional embedding scheme, the partition of electrons is an extra degree of freedom that cannot be automatically optimized but has to be determined a priori. In this work, as we are trying to reproduce the electronic structure of the central atom, we will use the Bader charges⁴⁶ of the central atom as a gauge to determine the partition. Basically, we optimize the embedding potentials with different electron partitions and conduct embedded cluster calculations using a consistent level of theory (PBE-PW). Then, we obtain the Bader charge of the central atom from the embedded cluster calculation and compare it with the supercell value. We will show that determining the correct charge for the cluster is critical in this study, and that the Bader charges indeed serve as a good gauge for this purpose.

Through the analysis described above, we find that the optimal cluster configurations are [SnZn₁₂S₄]¹⁸⁻ and [CuZn₁₂S₄]⁰, respectively (see Table I). In order to confirm that the embedded cluster models are able to reproduce the correct physics in the bulk, we conduct embedded cluster calculations of defect states with a consistent level of theory (PBE-PW) and compare the projected density of states (PDOSs) at the central atom with results from non-embedded bulk crystal supercell calculations of the same defect (see Fig. 7 and Table I).