Van Der Waals Interactions in Ionic and Semiconductor Solids

van der Waals (vdW) energy corrected density-functional theory [Phys. Rev. Lett. 102, 073005 (2009)] is applied to study the cohesive properties of ionic and semiconductor solids (C, Si, Ge, GaAs, NaCl, and MgO). The required polarizability and dispersion coefficients are calculated using the dielectric function obtained from time-dependent density-functional theory. Coefficients for ''atoms in the solid'' are then calculated from the Hirshfeld partitioning of the electron density. It is shown that the Clausius-Mossotti equation that relates the polarizability and the dielectric function is accurate even for covalently-bonded semiconductors. We find an overall improvement in the cohesive properties of Si, Ge, GaAs, NaCl, and MgO, when vdW interactions are included on top of the Perdew-Burke-Ernzerhof or Heyd-Scuseria-Ernzerhof functionals. The relevance of our findings for other solids is discussed. Cohesion in ionic and semiconductor solids arises mainly from electrostatic interactions and covalent bonding [1], and it is typically assumed that van der Waals (vdW) interactions play a minor role. Nevertheless, the qualitative and quantitative role of vdW interactions in solids has been a matter of discussion for quite some time [2–8]. Density-functional theory (DFT) is the method of choice for studying the bonding properties of solids. However, the widely used local-density approximation (LDA), generalized gradient approximation (GGA), as well as (screened) hybrid functionals are lacking the long-range vdW energy tail. Direct ab initio estimates of the vdW energy in solids have so far proven challenging. Recent many-body calculations using the random-phase approximation (RPA) [8–10], which includes the vdW energy seamlessly and accurately, yield significantly improved cohesive properties for a wide variety of solid-state systems over state-of-the-art (semi)-local DFT approximations. However, the vdW energy can only be rigorously defined in the large-distance limit. The aim of this Letter is to calculate accurate atomic polarizabilities, vdW coefficients , and vdW radii in ionic and semiconductor solids and to estimate the contribution of the long-range vdW energy to the cohesive properties. We find an overall improvement in the description of cohesion for model ionic and semiconductor solids when vdW interactions are included on top of the Perdew-Burke-Ernzerhof (PBE) [11] and the Heyd-Scuseria-Ernzerhof (HSE) [12,13] functionals. The first step to calculate the vdW energy in a solid is to determine the frequency-dependent polarizability per crystal unit cell. Two complementary schemes are followed in this work: (i) compute the polarizability per volume from the dielectric function of the solid, (ii) compute the polarizability per atom from cluster …

Cohesion in ionic and semiconductor solids arises mainly from electrostatic interactions and covalent bonding [1], and it is typically assumed that van der Waals (vdW) interactions play a minor role. Nevertheless, the qualitative and quantitative role of vdW interactions in solids has been a matter of discussion for quite some time [2][3][4][5][6][7][8]. Density-functional theory (DFT) is the method of choice for studying the bonding properties of solids. However, the widely used local-density approximation (LDA), generalized gradient approximation (GGA), as well as (screened) hybrid functionals are lacking the long-range vdW energy tail. Direct ab initio estimates of the vdW energy in solids have so far proven challenging. Recent many-body calculations using the random-phase approximation (RPA) [8][9][10], which includes the vdW energy seamlessly and accurately, yield significantly improved cohesive properties for a wide variety of solid-state systems over state-of-the-art (semi)-local DFT approximations. However, the vdW energy can only be rigorously defined in the large-distance limit. The aim of this Letter is to calculate accurate atomic polarizabilities, vdW coefficients, and vdW radii in ionic and semiconductor solids and to estimate the contribution of the long-range vdW energy to the cohesive properties. We find an overall improvement in the description of cohesion for model ionic and semiconductor solids when vdW interactions are included on top of the Perdew-Burke-Ernzerhof (PBE) [11] and the Heyd-Scuseria-Ernzerhof (HSE) [12,13] functionals.
The first step to calculate the vdW energy in a solid is to determine the frequency-dependent polarizability per crystal unit cell. Two complementary schemes are followed in this work: (i) compute the polarizability per volume from the dielectric function of the solid, (ii) compute the polarizability per atom from cluster extrapolation. We find good agreement between these two different approaches, which justifies the concept of ''atom in the solid.'' The atomic polarizability in binary (or more complex) crystals is obtained from the Hirshfeld partitioning of the electron density in a solid.
The frequency-dependent polarizability is related to the dielectric function as [14] where ð!Þ is the polarizability per volume V at frequency !, "ð!Þ is the dielectric function of the solid, and L is the so-called Lorentz factor, which depends on the ionicity and hybridization [14]. The constant L is known for two limits: (i) electron gas, where L ¼ 0, (ii) ideal rare-gas solids and ideal ionic crystals, where L ¼ 4=3. In the latter case, Eq. (1) is known as the Clausius-Mossotti (CM) equation.
In rare-gas and purely ionic crystals, the CM equation is an accurate approximation. However, the hybridization between nearest neighbors in a covalently-bonded semiconductor is thought to lead to a significantly reduced value of the Lorentz factor [15]. An alternative and direct way to compute the polarizability in solids is through a cluster extrapolation (CE). For ionic crystals, the CE and CM procedures yield essentially the same result [16]. Thus, we concentrate here on showing that CE and CM also lead to very similar polarizabilities for semiconductor crystals. We illustrate our case for Si, however, the same conclusions can be made for C and Ge [17]. In the case of hydrogen-saturated semiconductor clusters, denoted as X i H j (X ¼ C, Si, or Ge; see Fig. 1), there are four types of fourfold coordinated X atoms, the ''bulklike'' X atom, and surface X atoms bonded to one, Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. two, or 3 H atoms, respectively. For the X i H j cluster, the total static polarizability is partitioned as where n i and n j are the numbers of X and H atoms of the cluster, and X and H denote average static polarizabilities. To obtain a converged bulk silicon atom polarizability Si and a dispersion coefficient C SiSi 6 from the CE, we use a set of 30 Si i H j clusters, from SiH 4 to Si 172 H 120 , with the geometries taken from Ref. [18]. We calculated the frequency-dependent polarizability for every cluster at the time-dependent density-functional level of theory (TDLDA) using the same computational setup as in Ref. [18]. To test the accuracy of the additivity approximation in Eq. (2), we perform a least-squares fit of Si and H to a set of N ( ! 2) clusters, as shown in Fig. 1. For relatively small clusters, there are simply not enough bulk Si atoms, thus the obtained polarizabilities do not represent the crystal environment. However, for N ! 10 the fit converges to Si ¼ 26:8 AE 0:5 bohr 3 . Essentially the same converged value ( Si ¼ 26:6 AE 0:3 bohr 3 ) is obtained when starting from the largest cluster. These results indicate that the applied fitting procedure is robust. To obtain the bulk C SiSi 6 vdW coefficient, exactly the same procedure was applied. In Fig. 1 we also show the results of fitting C SiSi 6 on different clusters. The behavior of the fit is similar to the case of static polarizability. We note that a similar procedure was recently used by Grimme et al. to determine atomic C 6 coefficients from molecular TDDFT calculations [19].
In order to compare the CM and CE approaches, we carried out periodic and cluster TDDFT calculations for diamond, Si, and Ge as a function of the unit cell volume.
The computational details are reported in the supplemental material [17]. The bulk dielectric function was computed using both TDLDA and TDHSE with the so-called Nanoquanta (NQ) nonlocal exchange-correlation kernel [20], using the same approach as in Ref. [21]. When combined with TDHSE, the NQ kernel can correctly capture excitonic effects, which are known to play an important role for the optical spectrum of semiconductors [22]. TDHSE þ NQ method yields excellent agreement for the static dielectric constant " 1 with the experimental reference for ionic and semiconductor solids (see Table I and Ref. [21]).
The comparison of C SiSi 6 of CM and CE is shown in Fig. 2 as a function of the lattice constant. The static polarizability follows the same trend as C 6 . The corresponding values at the equilibrium volume of C, Si, and Ge are given in the supplemental material [17]. For the range of investigated unit cell volumes, Si and C SiSi 6 exhibit close to linear dependence on the lattice parameter. The maximum deviation between CM and CE is less than 0:5 bohr 3 for Si (2%), and 9:0 hartree bohr 6 for C SiSi 6 (6%). Considering that both CM and CE rely on a number of approximations, the agreement between them is striking. It is noteworthy that the polarizability and C 6 coefficients in semiconductor crystals are significantly reduced compared to the free-atom values (C SiSi 6;free ¼ 305 hartree bohr 6 ): at the equilibrium volume, the C 6 coefficient is reduced by roughly a factor of 2. In fact, this finding holds for all solids studied in this work, as can be seen in Table I by comparing the C 6 coefficients per unit cell with the respective freeatom sum. Seeking to understand the success of the CM equation for a wide variety of solids, beyond ''simple'' rare-gas and ideal ionic crystals, we use the DFT þ vdW method [24] to calculate the C 6 coefficients for atoms in semiconductors. The DFT þ vdW method accurately takes into account the short-range screening and hybridization contributions to the polarizability. However, the long-range screening is TABLE I. Calculated static dielectric constant " 1 and C 6 dispersion coefficients (in hartree bohr 6 per unit cell) at the experimental lattice constants using the TDHSE þ NQ method (fourth column-computed from the experimental dielectric function). The experimental dielectric constants are taken from Ref. [23]. The sixth column shows the sum of free-atom C 6 coefficients per unit cell. The last two columns show partitioned atomic C 6 coefficients used for the vdW energy correction. ¼ 280 hartree bohr 6less than 10% decrease from the corresponding free-atom values. Thus, we conclude that the main factor for the polarizability reduction in solids is the long-range electrostatic screening. In fact, Eq. (1) relates the dielectric function to the atomic polarizability by including the collective electrostatic screening from all the dipoles in the crystal. For crystals with cubic symmetry (which holds for all solids studied in this work), this results in L ¼ 4=3, yielding the CM equation. The constant L is only modified due to short-range screening effects. Since we found the short-range screening to play a minor role for C, Si, Ge, and GaAs, this explains the success of the CM equation for semiconductors. The inclusion of screening effects beyond hybridization in the DFT þ vdW method is the subject of our current work [25].
The good agreement between the CM and CE approaches allows us to use either C 6 ðVÞ values for calculating the vdW energy in solids. However, the best agreement with the experimental reference is achieved with TDHSE þ NQ calculations. In fact, when using the experimental dielectric function of Si, the C SiSi 6 coefficient agrees with the TDHSE þ NQ method to 3%. In the case of Ge, the TDHSE þ NQ dielectric constant " 1 is in excellent agreement with measurements of Cardona et al. [26], thus we use the polarizabilities and C 6 coefficients from the TDHSE þ NQ scheme for further discussion.
The missing long-range vdW energy is added to the DFT total energy on top of PBE [11] and HSE06 [13] functionals by using the DFT þ vdW approach [24] where R AB is the distance between atoms A and B, C AB 6 is the corresponding vdW coefficient, and R 0 A and R 0 B are the vdW radii. The latter are defined as R solid vdW ¼ ð solid = free Þ 1=3 R free vdW , where R free vdW is calculated as described in Ref. [24]. The damping function f d cuts off the interaction at short covalent distances. We use a Fermi-type damping function [27], , and d and s R are free parameters. The s R parameter scales the vdW radius and is adjusted for every DFT functional [28]. The d parameter controls the steepness of the damping function. For PBE we have determined s R ¼ 0:94 and d ¼ 20 on a database of accurate quantum-chemical calculations. For the HSE06 functional, s R ¼ 0:96 and d ¼ 20. The cohesive properties of solids show little sensitivity to a reasonable variation of up to 5% for s R and 15% for d.
To extend the DFT þ vdW scheme to solids with several atomic elements per unit cell, we define , where V Hirsh is the Hirshfeld volume of an atom A or B in the solid [24]. Using the additivity ansatz for the polarizability, we obtain B ði!Þ ¼ solid ði!Þ=ð AB þ 1Þ and A ði!Þ ¼ solid ði!Þ AB =ð AB þ 1Þ, where solid ði!Þ comes from Eq. (1). With this definition the C 6 coefficients (C AA 6 , C BB 6 , C AB 6 ) can be obtained from the Casimir-Polder integral using A ði!Þ and B ði!Þ. This scheme can be easily extended to ternary and more complex compounds by defining AC and so on. We have used the TDHSE þ NQ dielectric function at the experimental lattice constant to compute solid ði!Þ for GaAs, NaCl, and MgO. Figure 3 summarizes the errors in the cohesive properties obtained by using the Birch-Murnaghan equation of state [29] with standard and vdW-corrected PBE and HSE06 functionals for all solids studied in this work. The obtained values for every solid are shown in Tables IV and V of the supplemental material [17]. The zero-point energy (ZPE) within the harmonic approximation is added to the ground-state energies at each volume for every functional. All DFT calculations, except HSE06, were done with the FHI-aims all-electron code [30], with LDA, PBE, and PBEsol cohesive properties in excellent agreement with WIEN2K [31]. The HSE06 calculations were done with the VASP code [32], with ZPE correction added using the PBE phonon spectra.
For the traditional functionals, the considered cohesive properties follow the well-known trend: the lattice constants with LDA are too small, the bulk moduli and cohesive energies are too large; PBE shows the opposite tendency, and PBEsol lies in between of LDA and PBE. The HSE06 method partially cures the self-interaction problem in LDA and GGA functionals, leading to significantly better electronic structure (e.g., band gaps and dielectric constants) of solids compared to LDA and GGA functionals [33]. For the same reason, the HSE06 method also yields improved results for the lattice constants and bulk moduli of semiconductors and ionic solids. However, lattice constants are still somewhat overestimated, while the bulk moduli are underestimated in HSE06 calculations. Furthermore, the HSE06 scheme leads to underestimated cohesive energies (see supplemental material [17]), typically very similar to the PBE functional. When the vdW energy is coupled with the PBE functional, the errors in all cohesive properties are reduced by a factor of 2 when compared to experimental data. The same improvement is found for the HSE06 þ vdW method, except for the bulk moduli, where the errors remain roughly the same (underestimation for HSE06 calculations, overestimation for HSE06 þ vdW calculations). It is well known that the optimal value of the range separation parameter in the HSE06 method depends on the electrostatic screening properties of the solid [13,33]. Hence, the improvement of bulk moduli in the HSE06 þ vdW method requires further careful investigation.
For diamond, both the PBE þ vdW and the HSE06 þ vdW approaches yield slight overbinding. It is known, however, that anharmonic zero-point energy is important for diamond due to its relatively small mass, and will reduce the overbinding [34]. For Si and Ge, the vdW energy contributes $8% to the cohesive energy for both the PBE þ vdW and HSE06 þ vdW schemes, thus it cannot be considered negligible. The vdW energy also makes a visible contribution to the cohesive properties of GaAs, NaCl, and MgO. Both PBE þ vdW and HSE06 þ vdW results are in consistently better agreement with experiment than the standard functionals. The vdW energy contributes around 0:2 eV=atom to the cohesive energy for GaAs, NaCl, and MgO, and 9-16 GPa to the bulk modulus.
The vdW concept is only cleanly defined at distances where the atomic electron densities do not overlap between two atoms. For all values of s R and d in the damping function [Eq. (4)] used in this work, the interaction between the nearest neighbors in C, Si, and Ge crystals vanishes. The second-neighbor interaction is also damped depending on the s R parameter. For s R ¼ 1:0, the cohesive energies and the bulk moduli change by less than a few percent, while the lattice constants are virtually unmodified. Thus, we conclude that the obtained improvements with the PBE þ vdW and HSE06 þ vdW methods are due to the long-range vdW energy, and are not an artifact of the chosen damping function.
In summary, we obtained accurate vdW coefficients for atoms in ionic and semiconductor crystals using TDDFT calculations and the Clausius-Mossotti relation between the polarizability and the dielectric function. The DFT þ vdW method leads to overall improvement of the cohesive properties of ionic (NaCl, MgO) and semiconductor (Si, Ge, GaAs) solids. Despite the good agreement of the HSE06 þ vdW method with experiment, there still remains room for improvement. For example, nonadditive many-body vdW energy is not accounted for in the HSE06 þ vdW method. Also, anharmonic zero-point and temperature contributions to the cohesive properties should be included for the final comparison between experiment and theory. Note that it is well known that most GGA and hybrid functionals consistently underestimate cohesive energies and bulk moduli, and overestimate lattice constants for a wide variety of semiconductors, ionic solids, and metals [35]. We thus conclude that our findings about the importance of the long-range vdW energy are likely to be valid beyond the benchmark semiconductors and ionic solids studied in this work.
A. T. acknowledges support from the European Research Council (ERC Starting Grant VDW-CMAT).