{'fraction_non_alphanumeric': 0.03771073646850044, 'fraction_numerical': 0.013182912916719483, 'mean_word_length': 4.8311160384331115, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dielectric screening in extended systems using the self-consistent Sternheimer equation and localized basis sets

We develop a first-principles computational method for investigating the dielectric screening in extended systems using the self-consistent Sternheimer equation and localized non-orthogonal basis sets. Our approach does not require the explicit calculation of unoccupied electronic states, only uses two-center integrals, and has a theoretical scaling of order O(N^3). We demonstrate this method by comparing our calculations for silicon, germanium, diamond, and LiCl with reference planewaves calculations. We show that accuracy comparable to planewaves calculations can be achieved via a systematic optimization of the basis set.

Quasiparticle calculations based on the GW method [1][2][3][4][5][6] and calculations of optical spectra using the Bethe-Salpeter equation [6][7][8] have emerged as state-of-theart first-principles computational methods for the study of electronic and optical excitations both in extended solids and in nanoscale systems. While in many cases these methods can reproduce and even predict experimental photoemission and optical data with remarkable accuracy, the associated computational complexity limits their application to relatively small systems, typically in the range of ten to a hundred atoms. The main bottleneck in these calculations is the evaluation of the irreducible polarization propagator in the random-phase approximation (RPA), which is used for constructing the screened Coulomb interaction of the GW self-energy, and the electron-hole interaction kernel in the Bethe-Salpeter equation. The calculation of this propagator is performed using an expansion over unoccupied Kohn-Sham states [3][4][5][6][9][10][11], and has a theoretical scaling of O(N 4 ), N being the number of atoms in the system [12] [13]. In addition the polarization propagator is commonly calculated by expanding the Kohn-Sham states using planewaves basis sets [14]. While planewaves have several obvious advantages, large systems can become intractable due to the very large basis size. These observations clearly point to the need for new methods to study the dielectric screening in complex systems which (i) do not require the explicit calculation of unoccupied electronic states, (ii) use small basis sets, and (iii) exhibit a more favorable scaling with system size.
In order to reduce the computational complexity of Sternheimer-type approaches it appears advantageous to adapt these methods to the case of localized basis sets. Indeed, electronic structure codes based on localized representations are successfully employed to study very large systems, up to several thousands of atoms [29][30][31]. The key advantage of these basis sets is that they exploit the electron localization in the insulating state [32] in order to obtain a sparse representation of the ground-state density matrix. However, since these basis sets are optimized for providing an accurate description of the occupied Kohn-Sham manifold, it is not clear a priori what their performance would be in the case of excited state calculations.
In this work we demonstrate a first-principles pseudopotential scheme for calculating the inverse dielectric matrix in extended systems using non-orthogonal localized basis sets. Our scheme includes both localfield effects and the frequency-dependence of the inverse dielectric matrix, and does not require the explicit calculation of unoccupied Kohn-Sham states. This is achieved by adapting the self-consistent Sternheimer method of Ref. [12] to the case of a basis set of localized pseudo-atomic orbitals as implemented in the SIESTA code [29,33,34]. Here we illustrate the formalism and report calculations of the inverse dielectric matrix for prototypical semiconductors and insulators. We perform a systematic comparison with reference planewaves calculations, and we investigate the convergence of our results with the size of the local orbital basis set. A detailed report with extensive benchmarks and implementation details is provided elsewhere [35].
We calculate the inverse dielectric matrix in the random-phase approximation following Ref. [12]. The frequency-dependent variation of the density matrix ∆n(r, r , ω) corresponding to the non-local bare Coulomb potential v(r, r ) is given by: where ψ v is an occupied Kohn-Sham state of energy v and spin-unpolarized systems are considered for simplicity. The variations ∆ψ v (r, r , ±ω) of the single-particle states are solutions of the Sternheimer equations: where Kohn-Sham HamiltonianĤ and the projector P v on the occupied manifold act on the variable r , and W (r, r , ω) is the RPA frequency-dependent screened Coulomb interaction. In order to obtain the screened Coulomb interaction we use the change ∆V H of the Hartree potential resulting from the variation of the density matrix: Equations (1)-(4) need to be solved self-consistently. In the following we use this formalism in order to directly calculate the inverse dielectric matrix −1 (r, r , ω). This is done by replacing W in Eq. (2) by −1 , and by using the following relation instead of Eq. (4): In this case the self-consistent calculation is started by initializing the inverse dielectric matrix to the Dirac delta δ(r, r ). A derivation of the connection with the sumover-states approach and details on this formalism can be found in Refs. [12,35] respectively. We now move to a localized basis representation. We first make the choice of treating the variables r and ω in the Sternheimer equation Eq. (2) as parameters. The variable r will be represented on a real-space grid. Then we expand the quantities dependent on r in the basis of local orbitals φ i (r ): By replacing Eqs. (6), (7) in Eq. (2) and projecting both sides onto φ j (r ) we obtain the Sternheimer equation in matrix representation: Here H, S and ρ are the Kohn-Sham Hamiltonian, the overlap, and the density matrices, respectively. The matrix elements of ∆V [r,ω] in Eq. (8) are defined as: The solution of Eq. (8) yields the coefficients ∆c ± v [r,ω] which are used to construct the variation of the density matrix from Eq. (1): At this point it is possible to explicitly calculate ∆n(r, r , ω) on a real-space grid, and then perform the integral in Eq. (3) on the same grid. Equation (4) is also evaluated on the real space grid, and the resulting potential is projected on the local orbital basis using Eq. (9). This procedure is repeated until self-consistency is achieved, and is carried out independently for each value of the parameters r and ω.
For the purpose of comparison with standard planewaves calculations it is convenient to specialize this formalism to the case of crystalline solids. We proceed as follows: we incorporate the Bloch phase factors in the basis orbitals φ ik (r ) where k is a wavevector in the Brillouin zone, and R is a lattice vector. The orbitals φ i in Eq. (11) are the same as in Eq. (6), except that now the index i runs over the orbitals spanning one unit cell only. Equation (11) defines a basis of periodic functions and is used to expand the cell-periodic part u vk of the Bloch states ψ vk (r ) = e ik·r u vk (r ). By using a similar expansion for the variations of the wavefunctions in Eq. (7), and by projecting both sides of Eq. (2) on the basis orbitals, we obtain the Sternheimer equation for periodic systems: which is analogous to Eq. (8). In this case the Hamiltonian, overlap, and density matrices, as well as the coefficient vectors, are all resolved in momentum space. A detailed derivation of these equations is provided in Ref. [35].
In the present formalism the dependence on the realspace variable r is represented on the local orbital basis, while the dependence on r is represented on a realspace grid. This choice carries the following advantages: (i) we do not use a product-basis expansion for non-local quantities, therefore we avoid issues related to the representability of the screened Coulomb interaction in localized basis sets. (ii) We only have two-center integrals in our formulation, therefore we do not need the three and four-center integrals arising in product basis expansions [10,11,[36][37][38]]. (iii) The self-consistent calculation of the inverse dielectric matrix avoids from the outset the issues associated with the inversion of response functions represented in local orbital basis sets [39]. (iv) Since Eq. (8) can be solved using sparse linear algebra and needs to be performed for every occupied state and every r on the real-space grid, this method has a theoretical scaling of O(N 3 ).
The key approximation in our formulation is the expansion of the variation ∆ψ ± v[r,ω] of the single-particle states in the basis of local orbital [cf. Eq. (7)]. In fact local orbital basis sets are typically optimized to accurately describe the occupied states manifold, while the variation of the density matrix arises from the components of ∆ψ ± v[r,ω] in the manifold of unoccupied states [25]. This clearly points to the need of carefully optimizing the local orbital basis sets for calculations of the dielectric screening. A systematic assessment of the performance of multiple-ζ polarized pseudo-atomic orbital basis sets is provided in Ref. [35].
Our implementation is based on the SIESTA code, and uses a basis of strictly localized numerical pseudo-atomic orbitals [29,33,34]. In the following we discuss benchmark results for silicon, diamond, germanium and LiCl. Calculations are performed within the local-density approximation (LDA) to density-functional theory [40,41], and with norm-conserving pseudopotentials [42]. The inverse dielectric matrices are calculated by sampling the Brillouin zone using a shifted 10×10×10 grid for diamond, silicon and LiCl, and a shifted 12×12×12 grid for germanium. We use the lattice parameters 5.43Å, 3.56 A, 5.65Å, and 5.13Å for silicon, diamond, germanium, and LiCl, respectively [43]. Using a standard triple-ζ polarized (TZP) basis we obtain direct band gaps of 2.55 eV (Si), 5.60 eV (diamond), 0.04 eV (Ge), 5.99 eV (LiCl), in line with standard LDA results. In order to generate basis sets for silicon with a large number of ζ's we use an energy-shift parameter [34] of 10 meV. This shift leads to localization radius of 9.3Å for silicon which is slightly larger than those adopted in standard ground-state calculations using SIESTA. We compare our results with planewaves calculations performed using ABINIT [44] and YAMBO [45] software packages, with same the pseudopotentials and Brillouin-zone sampling for consistency [46]. Planewaves calculations are carried out using kinetic energy cutoffs of 20 Ry for Si and Ge, and of 60 Ry for diamond and LiCl. The RPA dielectric matrices are calculated using 92 conduction bands and kinetic energy cutoffs of 6.9 Ry for silicon, germanium and LiCl, and of 12 Ry for diamond. Figure 1 shows the calculated static macroscopic dielectric constant of silicon as a function of basis set size. The size of the basis is increased by including additional ζ's using the split-norm procedure, as well as polarization orbitals [34]. We find that the single-ζ polarized (SZP) basis already gives results which are within 10% of the reference planewaves calculation. This suggests that a suitably optimized basis with a few orbitals should be able to yield converged results. We also observe a mono- tonic convergence with the number of basis orbitals, and a consistently superior performance of the basis sets including polarization orbitals (Fig. 1).
We can rationalize the convergence of the dielectric constant with basis size shown in Fig. 1 by considering the Penn model for the dielectric constants of semiconductors [47]. In this simple model the dielectric constant is related to the Penn gap E P and to the plasma frequency ω p by 0 = 1 + ( ω p /E P ) 2 . In our calculations we extract the Penn gap for each basis set considered in Fig. 1 using the main peak in the joint density of states (JDOS) [48] [49]. Figure 2 shows that 0 − 1 is approximately a linear function of 1/E 2 P . In particular the values calculated using the unpolarized basis sets nicely fall on the line obtained using the plasma frequency of silicon ω p =16.6 eV. For comparison, the inset of Fig. 2 shows that the calculated dielectric constants do not correlate with the direct band gap of silicon. These observations lead us to conclude that the Penn gap, and more generally the low-energy structure of the JDOS, is a good indicator of the quality of the basis set for calculating dielectric constants. This is helpful for a quick assessment of the quality of a basis set without explicitly calculating the screening. Table I shows the first few elements of the inverse dielectric matrix −1 GG (q → 0, ω = 0) for several local orbital basis sets. We observe that the calculated values for the first few shells compare well with the reference planewaves calculation, similarly to the head of the dielectric matrix. For higher shells, however, the relative agreement of the basis sets worsens slightly. This may indicate a systematic deficiency in the local orbital ba- sis in case of the higher Fourier components of dielectric matrix, and clearly deserves further investigation. Figure 3 shows the frequency-dependent dielectric function of silicon calculated using our method, together with a reference planewaves calculation. We observe that, as expected, the performance of the minimal singleζ basis is poor. In fact, spectral weight is incorrectly transferred from the main absorption peak to higher energies. On the contrary, the TZP basis (17 orbitals per atom) performs quite well, with all the main features of the planewaves spectra correctly reproduced. We still observe however some small transfer of spectral weight and a slight blueshift of the high-energy peaks.
In order to demonstrate the generality of of our approach we present in Tab. II our calculated dielectric constants for germanium, diamond, and LiCl. The values shown in Tab. II are obtained by using standard SIESTA   II. Dielectric constants of Si, Ge, diamond, and LiCl calculated using our Sternheimer approach, and compared to reference planewaves calculations and experimental data. We report both our results obtained using the minimal SZ basis and the TZP basis, as generated using the standard settings of SIESTA.
settings for the basis [50]. Table II provides further support to our previous finding by showing that the TZP basis provides results which lie within 1-5% of the reference planewaves calculations. We expect further improvement upon designing basis sets specifically optimized for the Sternheimer scheme proposed in this work. In particular the use of numerical diffuse orbitals [51,52] deserves a systematic assessment.
In conclusion, we have introduced and demonstrated a method for calculating the inverse dielectric matrix of ex-tended systems combining localized non-orthogonal basis sets with the self-consistent Sternheimer equation. Our method does not require the calculation of unoccupied electronic states, does not require the explicit inversion of the dielectric matrix, uses only two-center integrals, and has a theoretical scaling with system size of O(N 3 ). Our implementation based on the pseudo-atomic orbitals of the SIESTA code shows that results with accuracy comparable to planewaves calculations can be obtained upon optimization of the localized basis set. Our formulation is completely general and there should be no difficulties in adapting our method to other local basis implementations. We believe this work represents a stepping stone towards quasiparticle GW calculations for large and complex systems using localized basis sets.