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 . 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.
pacs:71.15.-m 71.15.Ap 71.15.Dx
Quasiparticle calculations based on the method Hedin (1965); Hedin and Lundqvist (1969); Hybertsen and Louie (1986); Aryasetiawan and Gunnarsson (1998); Aulbur et al. (2000); Onida et al. (2002) and calculations of optical spectra using the Bethe-Salpeter equation Onida et al. (1995); Rohlfing and Louie (1998); Onida et al. (2002) have emerged as state-of-the-art 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 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 Hybertsen and Louie (1986); Rohlfing et al. (1993, 1995); Aryasetiawan and Gunnarsson (1998); Aulbur et al. (2000); Onida et al. (2002); Blase et al. (2011), and has a theoretical scaling of , being the number of atoms in the system Giustino et al. (2010) 111By using planewaves basis sets and fast Fourier transform techniques a scaling can be obtained.. In addition the polarization propagator is commonly calculated by expanding the Kohn-Sham states using planewaves basis sets Ihm et al. (1979). 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.
Several methods have been proposed in order to avoid the explicit calculation of unoccupied electronic states in the study of electronic excitations Levine and Allan (1989, 1991); Reining et al. (1997); Andrade et al. (2007); Wilson et al. (2008); Bruneval and Gonze (2008); Wilson et al. (2009); Rocca et al. (2010); Giustino et al. (2010); Umari et al. (2010); Berger et al. (2010). Some of these methods rely on the Sternheimer equation Levine and Allan (1989, 1991); Reining et al. (1997); Andrade et al. (2007); Giustino et al. (2010), which is routinely used in density-functional perturbation theory Baroni et al. (2001), and closely related to the coupled-perturbed Hartree-Fock method found in the quantum chemistry literature McWeeny (1962); Gerratt (1968); Liang et al. (2005). In particular, in Ref. Giustino et al. (2010) the Sternheimer method is used self-consistently in order to calculate the full frequency-dependent inverse dielectric function in extended systems, and an empirical pseudopotential implementation is described as a proof-of-concept.
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 Soler et al. (2002); Skylaris et al. (2005); Bowler and Miyazaki (2010). The key advantage of these basis sets is that they exploit the electron localization in the insulating state Kohn (1996) 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 local-field 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. Giustino et al. (2010) to the case of a basis set of localized pseudo-atomic orbitals as implemented in the SIESTA code Sánchez-Portal et al. (1996); Artacho et al. (1999); Soler et al. (2002). 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 Hübener et al. (2012).
We calculate the inverse dielectric matrix in the random-phase approximation following Ref. Giustino et al. (2010). The frequency-dependent variation of the density matrix corresponding to the non-local bare Coulomb potential is given by:
where is an occupied Kohn-Sham state of energy and spin-unpolarized systems are considered for simplicity. The variations of the single-particle states are solutions of the Sternheimer equations:
where Kohn-Sham Hamiltonian and the projector on the occupied manifold act on the variable , and is the RPA frequency-dependent screened Coulomb interaction. In order to obtain the screened Coulomb interaction we use the change 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 . This is done by replacing in Eq. (2) by , 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 . A derivation of the connection with the sum-over-states approach and details on this formalism can be found in Refs. Giustino et al. (2010); Hübener et al. (2012) respectively.
We now move to a localized basis representation. We first make the choice of treating the variables and in the Sternheimer equation Eq. (2) as parameters. The variable will be represented on a real-space grid. Then we expand the quantities dependent on in the basis of local orbitals :
Here , and are the Kohn-Sham Hamiltonian, the overlap, and the density matrices, respectively. The matrix elements of in Eq. (8) are defined as:
At this point it is possible to explicitly calculate 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 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
where is a wavevector in the Brillouin zone, and is a lattice vector. The orbitals in Eq. (11) are the same as in Eq. (6), except that now the index 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 of the Bloch states . 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. Hübener et al. (2012).
In the present formalism the dependence on the real-space variable is represented on the local orbital basis, while the dependence on is represented on a real-space 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 Rohlfing et al. (1995); Blase and Ordejón (2004); Koval et al. (2010); Foerster et al. (2011); Blase et al. (2011). (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 Gidopoulos and Lathiotakis (2011). (iv) Since Eq. (8) can be solved using sparse linear algebra and needs to be performed for every occupied state and every on the real-space grid, this method has a theoretical scaling of .
The key approximation in our formulation is the expansion of the variation 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 in the manifold of unoccupied states Baroni et al. (2001). 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. Hübener et al. (2012).
Our implementation is based on the SIESTA code, and uses a basis of strictly localized numerical pseudo-atomic orbitals Sánchez-Portal et al. (1996); Artacho et al. (1999); Soler et al. (2002). 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 Ceperley and Alder (1980); Perdew and Zunger (1981), and with norm-conserving pseudopotentials Troullier and Martins (1991). The inverse dielectric matrices are calculated by sampling the Brillouin zone using a shifted 101010 grid for diamond, silicon and LiCl, and a shifted 121212 grid for germanium. We use the lattice parameters 5.43 Å, 3.56 Å, 5.65 Å, and 5.13 Å for silicon, diamond, germanium, and LiCl, respectively Hybertsen and Louie (1987). 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 Artacho et al. (1999) 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 Gonze et al. (2009) and YAMBO Marini et al. (2009) software packages, with same the pseudopotentials and Brillouin-zone sampling for consistency 222 We note that there is a small difference between SIESTA and ABINIT in how the local part of the pseudopotential is constructed.. 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 Artacho et al. (1999). 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 monotonic 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 Penn (1962). In this simple model the dielectric constant is related to the Penn gap and to the plasma frequency by . 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) Cardona and Yu (1999) 333 For a reference basis set we choose the Penn gap as the energy of the main peak of the JDOS. The cumulative integral of the JDOS from the absorption threshold to this energy is then used to define unambiguously the Penn gap for the other basis sets.. Figure 2 shows that is approximately a linear function of . In particular the values calculated using the unpolarized basis sets nicely fall on the line obtained using the plasma frequency of silicon =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 1 shows the first few elements of the inverse dielectric matrix 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 basis 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. 2 our calculated dielectric constants for germanium, diamond, and LiCl. The values shown in Tab. 2 are obtained by using standard SIESTA settings for the basis 444In this table we use a split-norm of 0.15 and an energy shift of 200 meV for silicon, diamond and germanium, and end nergy shift of 70 meV for LiCl. This choice leads to small localization radii, for instance we obtain a radius of 6.4 Å in the case of Si.. Table 2 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 Ciofini and Adamo (2007); García-Gil et al. (2009) deserves a systematic assessment.
Ref. Kittel (1986).
Ref. Ashcroft and Mermin (1976).
In conclusion, we have introduced and demonstrated a method for calculating the inverse dielectric matrix of extended 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 . 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 calculations for large and complex systems using localized basis sets.
The authors would like to thank E. Artacho for fruitful discussions. This work is funded by the European Research Council under the European Community’s Seventh Framework Programme, Grant No. 239578, and Spanish MICINN Grant FIS2009-12721-C04-01.
- Hedin (1965) L. Hedin, Phys. Rev. 139, A796 (1965).
- Hedin and Lundqvist (1969) L. Hedin and S. Lundqvist, Solid State Physics, edited by F. Seitz, D. Turnbull, and A. Ehrenreich, Vol. 23 (Academic Press, New York, 1969) p. 1.
- Hybertsen and Louie (1986) M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
- Aryasetiawan and Gunnarsson (1998) F. Aryasetiawan and O. Gunnarsson, Rep. Prog. Phys. 61, 237 (1998).
- Aulbur et al. (2000) W. G. Aulbur, L. Jönsson, and J. W. Wilkins, Solid State Physics, edited by A. Ehrenreich and F. Spaepen, Vol. 54 (Academic, San Diego, 2000) p. 1.
- Onida et al. (2002) G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
- Onida et al. (1995) G. Onida, L. Reining, R. W. Godby, R. Del Sole, and W. Andreoni, Phys. Rev. Lett. 75, 818 (1995).
- Rohlfing and Louie (1998) M. Rohlfing and S. G. Louie, Phys. Rev. Lett. 81, 2312 (1998).
- Rohlfing et al. (1993) M. Rohlfing, P. Krüger, and J. Pollmann, Phys. Rev. B 48, 17791 (1993).
- Rohlfing et al. (1995) M. Rohlfing, P. Krüger, and J. Pollmann, Phys. Rev. B 52, 1905 (1995).
- Blase et al. (2011) X. Blase, C. Attaccalite, and V. Olevano, Phys. Rev. B 83, 115103 (2011).
- Giustino et al. (2010) F. Giustino, M. L. Cohen, and S. G. Louie, Phys. Rev. B 81, 115105 (2010).
- (13) By using planewaves basis sets and fast Fourier transform techniques a scaling can be obtained.
- Ihm et al. (1979) J. Ihm, A. Zunger, and M. L. Cohen, J. Phys. C: Solid State Phys. 12, 4409 (1979).
- Levine and Allan (1989) Z. H. Levine and D. C. Allan, Phys. Rev. Lett. 63, 1719 (1989).
- Levine and Allan (1991) Z. H. Levine and D. C. Allan, Phys. Rev. B 43, 4187 (1991).
- Reining et al. (1997) L. Reining, G. Onida, and R. W. Godby, Phys. Rev. B 56, R4301 (1997).
- Andrade et al. (2007) X. Andrade, S. Botti, M. A. L. Marques, and A. Rubio, J. Chem. Phys. 126, 184106 (2007).
- Wilson et al. (2008) H. F. Wilson, F. Gygi, and G. Galli, Phys. Rev. B 78, 113303 (2008).
- Bruneval and Gonze (2008) F. Bruneval and X. Gonze, Phys. Rev. B 78, 085125 (2008).
- Wilson et al. (2009) H. F. Wilson, D. Lu, F. Gygi, and G. Galli, Phys. Rev. B 79, 245106 (2009).
- Rocca et al. (2010) D. Rocca, D. Lu, and G. Galli, J. Chem. Phys. 133, 164109 (2010).
- Umari et al. (2010) P. Umari, G. Stenuit, and S. Baroni, Phys. Rev. B 81, 115104 (2010).
- Berger et al. (2010) J. A. Berger, L. Reining, and F. Sottile, Phys. Rev. B 82, 041103 (2010).
- Baroni et al. (2001) S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Rev. Mod. Phys. 73, 515 (2001).
- McWeeny (1962) R. McWeeny, Phys. Rev. 126, 1028 (1962).
- Gerratt (1968) J. Gerratt, J. Chem. Phys. 49, 1719 (1968).
- Liang et al. (2005) W. Liang, Y. Zhao, and M. Head-Gordon, J. Chem. Phys. 123, 194106 (2005).
- Soler et al. (2002) J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón, and D. Sánchez-Portal, J. Phys. Condens. Matter 14, 2745 (2002).
- Skylaris et al. (2005) C.-K. Skylaris, P. D. Haynes, A. A. Mostofi, and M. C. Payne, J. Chem. Phys. 122, 084119 (2005).
- Bowler and Miyazaki (2010) D. R. Bowler and T. Miyazaki, J. Phys. Condens. Matter 22, 074207 (2010).
- Kohn (1996) W. Kohn, Phys. Rev. Lett. 76, 3168 (1996).
- Sánchez-Portal et al. (1996) D. Sánchez-Portal, E. Artacho, and J. M. Soler, J. Phys. Condens. Matter 8, 3859 (1996).
- Artacho et al. (1999) E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García, and J. M. Soler, Phys. Status Solidi B 215, 809 (1999).
- Hübener et al. (2012) H. Hübener, M. A. Pérez-Osorio, P. Ordejón, and F. Giustino, (2012), arXiv:1202.XXXX.
- Blase and Ordejón (2004) X. Blase and P. Ordejón, Phys. Rev. B 69, 085111 (2004).
- Koval et al. (2010) P. Koval, D. Foerster, and O. Coulaud, Phys. Status Solidi B 247, 1841â1848 (2010).
- Foerster et al. (2011) D. Foerster, P. Koval, and D. Sanchez-Portal, J. Chem. Phys. 135, 074105 (2011).
- Gidopoulos and Lathiotakis (2011) N. I. Gidopoulos and N. N. Lathiotakis, (2011), arXiv:1107.6007.
- Ceperley and Alder (1980) D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
- Perdew and Zunger (1981) J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
- Troullier and Martins (1991) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
- Hybertsen and Louie (1987) M. S. Hybertsen and S. G. Louie, Phys. Rev. B 35, 5585 (1987).
- Gonze et al. (2009) X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken, F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Caracas, M. Côté, T. Deutsch, L. Genovese, P. Ghosez, M. Giantomassi, S. Goedecker, D. Hamann, P. Hermet, F. Jollet, G. Jomard, S. Leroux, M. Mancini, S. Mazevet, M. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rignanese, D. Sangalli, R. Shaltaf, M. Torrent, M. Verstraete, G. Zerah, and J. Zwanziger, Comput. Phys. Commun. 180, 2582 (2009).
- Marini et al. (2009) A. Marini, C. Hogan, M. Grüning, and D. Varsano, Comput. Phys. Commun. 180, 1392 (2009).
- (46) We note that there is a small difference between SIESTA and ABINIT in how the local part of the pseudopotential is constructed.
- Penn (1962) D. R. Penn, Phys. Rev. 128, 2093 (1962).
- Cardona and Yu (1999) M. Cardona and P. Y. Yu, Fundamentals of Semiconductors: Physics and Materials Properties (Springer, 1999).
- (49) For a reference basis set we choose the Penn gap as the energy of the main peak of the JDOS. The cumulative integral of the JDOS from the absorption threshold to this energy is then used to define unambiguously the Penn gap for the other basis sets.
- (50) In this table we use a split-norm of 0.15 and an energy shift of 200 meV for silicon, diamond and germanium, and end nergy shift of 70 meV for LiCl. This choice leads to small localization radii, for instance we obtain a radius of 6.4 Å in the case of Si.
- Ciofini and Adamo (2007) I. Ciofini and C. Adamo, J. Phys. Chem. A 111, 5549 (2007).
- García-Gil et al. (2009) S. García-Gil, A. García, N. Lorente, and P. Ordejón, Phys. Rev. B 79, 075441 (2009).
- Kittel (1986) C. Kittel, Introduction to Solid State Physics (Wiley, New York, 1986) p. 207.
- Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Harcourt, New York, 1976) p. 553.