Hybrid functional study of structural, electronic, bonding and optical properties of CdSiP2

Computational Materials Science 117 (2016) 472–477

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Hybrid functional study of structural, electronic, bonding and optical properties of CdSiP2 Jianping Xiao a,b, Zhiyu He a,⇑, Shifu Zhu a, Baojun Chen a, Gang Jiang c a

Department of Material Science, Sichuan University, Chengdu 610064, China College of Electrical and Information Engineering, Southwest University for Nationalities, Chengdu 610041, China c Institute of Atomic and Molecular Physics, Sichuan University, Chengdu 610065, China b

a r t i c l e

i n f o

Article history: Received 26 July 2015 Received in revised form 28 January 2016 Accepted 8 February 2016 Available online 1 March 2016 Keywords: Hybrid functional Electronic structure Bonding properties Optical properties

a b s t r a c t The structural, electronic, bonding and optical properties of CdSiP2 have been investigated by hybrid density functional theory. The results of structural parameters and band gap calculated using HSE06 with 25% short-range HF exchange agree well with experimental data. The nature of the band gap in CdSiP2 is predicted to be pseudodirect by comparing the band structures of CdSiP2 and GaP. The bonding properties are analyzed based on the Bader charge and the electronic localization function. In addition, the optical spectrums are derived and main spectral features are interpreted in terms of the electronic structure. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction In lack of suitable all-solid-state laser sources, optical parametric down conversion employing a nonlinear optical (NLO) crystal is an available approach to cover the mid-infrared (mid-IR) spectral range for wavelength generation. At present, I–III–VI2 and II–IV–V2 chalcopyrite semiconductors, such as AgGaSe2, AgGaxIn1
xSe2, ZnGeP2, CdGeAs2, and CdSiP2, are the main IR NLO materials. Among these materials, CdSiP2 has recently attracted much more attention due to its large effective nonlinear coefficient [1], wide transparent region and good mechanical properties [2]. However, it is still a great challenge to prepare stoichiometric and defect-free CdSiP2 single crystals for the reason of the high melting point, extremely high vapor pressure, huge anisotropic thermal expansion and etc., during the growth process [3,4]. Up to now, due to lack of high-quality CdSiP2 crystals, the properties of this material have not been sufficiently investigated experimentally. In order to further understand fundamental properties, many theoretical investigations have been devoted to CdSiP2. Chiker et al. [5] performed the calculation using the FPLAPW method within LDA and reported that CdSiP2 has a pseudodirect band gap of 1.10 eV. DFT calculations (GGA) reported by Lv et al. [6] and He et al. [7] also underestimate the band gap. And their predictions ⇑ Corresponding author. Tel./fax: +86 28 85412745. E-mail address: [email protected] (Z. He). http://dx.doi.org/10.1016/j.commatsci.2016.02.014 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

about the nature of the gap are controversial. Shaposhnikov et al. [8] identified the gap as direct and gave reasonable gap values: 1.91 and 2.05 eV for GW0 approximation and modified Becke– Johnson (mBJ) exchange potential. We noted that the band gap values have a wide distribution and there is no consensus on the nature of the band gap. Thus, an accurate determination of the bands combined with zone folding arguments needs further discussion. Furthermore, considering that the bulk modulus is an important mechanical property for designing optical devices, the issue of divergences between the theoretical and experimental bulk modulus also should be addressed. All previous DFT calculations of the structural parameters were carried out within LDA or GGA. LDA always underestimates the equilibrium volume by 3–6% [9], while GGA overestimates the equilibrium volume by about 2–5% [10]. However, the bulk modulus is highly sensitive to the equilibrium volume [10]. Fortunately, the hybrid density functional which mixes a certain fraction of the exact Hartree–Fock (HF) exchange with the local or semi-local density functional has proved to be an efficient approach to remedy both the volume error and band gap error. As far as we know, there is no any theoretical study of CdSiP2 using hybrid density functional theory. Therefore, in this case, we have calculated the structural, electronic, bonding and optical properties using the HSE06 functional. Based on the obtained results, we have also discussed that whether the addition of 25% Hartree–Fock exchange in the hybrid density functional is appropriate.

473

J. Xiao et al. / Computational Materials Science 117 (2016) 472–477

2. Computational details All calculations were performed using the projector augmented wave (PAW) method, as implemented in the VASP package [11,12]. Both PBE [13] functional and hybrid density functional were employed as the exchange–correlation functionals. The valence electronic configurations for Cd, Si, and P were 4d105s2, 3s23p2, and 3s22p3, respectively. The Gaussian smearing of 0.02 eV was applied in all calculations except for the total energy and density of states (DOS) calculations. An energy cut-off of 400 eV for a plane-wave basis and a C-centered 4  4  4 k-point mesh were used for the structural optimization. The atomic relaxation was stopped when the force on each atom was less than 0.01 eV/Å. Total energy convergence tolerance was less than 1.0  10
5 eV/atom. DOS and dielectric properties were calculated with a denser k-point mesh of 6  6  6. In our hybrid density functional study, the HSE06 hybrid functional [14] was used. The exchange–correlation energy of HSE06 is given by following expression:

EHSE xc

¼a

lÞ þ ð1
a

EHF;SR ð x

ÞEPBE;SR ð x

lÞ þ

lÞ þ

EPBE;LR ð x

EPBE c

ð1Þ

The exchange interaction is separated into a short-ranged (SR) and long-ranged (LR) part. The slowly decaying LR part is totally represented by PBE exchange. The SR part is rational mixing of the HF and PBE exchange. For the HSE06 functional the screening parameter l which defines the range-separation is set to 0.2 Å
1. The mixing parameter a specifies the fraction of the HF exchange and PBE exchange. The default value a = 0.25 was used in our study. The VESTA package [15] was used to analyze the band decomposed charge densities and the electron localization.

3. Results and discussions 3.1. Structural parameters and bulk modulus CdSiP2 crystallizes in the body-centered tetragonal chalcopyrite  structure (space group I42d). The chalcopyrite structure can be treated as a superlattice by doubling the zinc-blende structure along the z-axis with small structural distortions. The structural parameters are the lattice constant a, tetragonal distortion parameter g = c/2a and anion displacement parameter u [16]. Here, u is determined by

u¼


2 1 1 c 1 2

2 32a2 16

ð2Þ

In the ideal structure, g = 1 and u = 1/4. The structure optimization was performed as follows: The shape and internal coordinates were optimized at each fixed volume. The total energy and equilibrium volume were obtained by fitting the energy–volume data to a third-order Birch–Murnaghan equation of state [17,18]. The total energy versus primitive cell volume calculated with PBE and HSE06 are presented in Fig. 1. The theoretical lattice constant a, tetragonal distortion parameter g, anion displacement parameter u, primitive cell volume (V0), bulk modulus (B0) and pressure derivative of bulk modulus ðB00 Þ together with available experimental values are summarized in Table 1. HSE06 reduces overestimation of the lattice constant as compared with PBE. The deviation decrease from 1% (PBE) to 0.3% (HSE06) for the parameter a. Tetragonal distortion parameter g and anion displacement parameter u are less sensitive to different functionals. As a consequence, HSE06 corrects the overestimation of equilibrium volume from 3% to 1% and the underestimation of bulk modulus roughly from 30% to 17%.

Fig. 1. Total energy versus primitive cell volume of CdSiP2 calculated within PBE (squares) and HSE06 (circles).

Table 1 Lattice constant a, tetragonal distortion parameter g, anion displacement parameter u, primitive cell volume V0, bulk modulus B0 and pressure derivative of bulk modulus B00 of CdSiP2 calculated using PBE and HSE06 compared to experimental values.

a b

Functional

a (Å)

g

u

V0 (Å3)

B0 (GPa)

B00

PBE HSE06 Expt.

5.738 5.699 5.679a

0.920 0.919 0.918a

0.2920 0.2924 0.2968a

173.75 170.14 168.26a

74 83 97b

4.51 4.36

Ref. [19]. Ref. [20].

3.2. Electronic properties The band structure of CdSiP2 calculated using HSE06 is shown in Fig. 2. The conduction band minimum (CBM) as well as the valence band maximum (VBM) is situated at the C point. In the ideal tetragonal structure, the primitive cell volume of the chalcopyrite (CH) structure is four times as large as that of the zincblende (ZB) structure, while the Brillouin zone of the chalcopyrite structure is a quarter of that one. As a consequence, the highsymmetry points X in the z direction, X in the x and y direction and L of the zinc-blende Brillouin zone fold back to C, T and N of the chalcopyrite Brillouin zone. Cubic GaP is identified as the corresponding III–V analogue of CdSiP2. In order to better understand the nature of the band gap of CdSiP2, band structures of GaP and CdSiP2 calculated within HSE06 are both plotted in Fig. 3. The eigenvalues at the

Fig. 2. Band structure of CdSiP2 along the high-symmetry line T–C–X–P–N–C calculated within HSE06. The zero of energy is set at the valence-band maximum.

474

J. Xiao et al. / Computational Materials Science 117 (2016) 472–477

states. And the hybridizations among s, p states of Si atoms and P atoms in a wide range can be observed. The d states of Cd atoms have little impact on the band structure near the Fermi level. The band decomposed charge density corresponding to the VBM and the CBM shown in Fig. 5 illustrates the origins of the VBM and CBM. As expected, the pz electrons distribute along the c-axis, whereas the (px, py) electrons distribute in the twodimensional ab-plane perpendicular to the c-axis. The charge density of the VBM (C4v) is dominated by the pz orbits of P atoms with a minor contribution from Cd atoms. As a matter of fact, the contribution from Cd atoms almost can be ignored. The charge density of the CBM (C3c) is derived from pz orbits of P atoms and (px, py) orbits of Si atoms. 3.3. Bonding properties Fig. 3. Band structures of CdSiP2 and its III–V analogue GaP calculated within HSE06. The zero of energy is set at the respective valence-band maximum.

high-symmetry points in the respective Brillouin zones are summarized in Table 2. For GaP, the CBM is situated at the X point and the X1c state is lower than the C1c state obviously. In the case of CdSiP2, the CBM is situated at the C point. The correspondence between CBMs of GaP and CdSiP2 indicates that the gap of CdSiP2 (C4v ? C3c) is derived from the indirect transition of GaP (C15v ? X1c). Based on this comparison, a conclusion was drawn that CdSiP2 has a so-called pseudodirect band gap. And it should be noted that the optical transition at C from the VBM to the CBM (C4v ? C3c) is forbidden by symmetry [21]. Due to reduce Td symmetry of the zinc-blende structure to D2d symmetry of the chalcopyrite structure, the triply degenerate C15v state splits into a doubly degenerate C5v and a non-degenerate C4v state in the Brillouin zone of CdSiP2. The crystal field splitting energy DCF is defined as the energy difference between the C4v and C5v state. DCF is negative when the C5v state lies below the C4v state. For CdSiP2, DCF =
0.27 eV follows the trend that most of chalcopyrite compounds have negative DCF when g is less than 1 [22]. The PBE functional predicts pseudodirect gap as 1.41 eV. On the other hand, the pseudodirect gap calculated using HSE06 is 2.06 eV which is much closer to the experimental results of 2.2–2.45 eV [23] than the PBE result. Furthermore, the band gap calculated by employing the highly accurate GW method is 2.10 eV, which also confirms our HSE result. The projected density of states (PDOS) of CdSiP2 calculated using HSE06 is depicted in Fig. 4. The lower valence bands (
14 eV 
11 eV) mainly consist of Cd-4d, Si-3s and P-2s states. The valence bands (
10 eV 
6 eV) are mainly made up of Cd-4d and Si-3s states. Cd-4d states at about
9 eV show a strong localization character. Near the Fermi level, the valence bands (
5 eV  0 eV) consist of Si-3p and P-2p states. Above the Fermi level, the conduction bands mostly derive from Si-3p and P-2p

We performed calculations of the electronic charge distribution with Bader charge analysis [26]. The results are listed in Table 3. Each Cd atom and Si atom loses about 0.59 electrons and 1.80 electrons to P atoms, respectively. The transfer charges are considerably smaller than the values of a pure ionic compound indicated that there is a covalent contribution to the bonding of CdSiP2. To examine bonding characters between different atoms, the electron localization function (ELF) contour plot is shown in Fig. 6. ELF values are restricted between 0 and 1. The upper limit ELF = 1 points to perfect localization, which can be interpreted as bonding pairs or lone pairs. ELF = 0.5 corresponds to electrongas-like pair probability [27,28]. The lower limit ELF = 0 corresponds to no localization. ELF is a quantitative measure to discriminate metallic, covalent and ionic bonding effectively. In Fig. 6, the high electron localization with ELF = 0.89 in the middle regions between the P atom and the Si atom may indicate significantly covalent character of the P–Si bonding. The high electron localization around the P atom and a narrow channel with a very low ELF value separating valence shells of the P atom and the neighboring Cd atom may indicate that P–Cd bonding is ionic. 3.4. Optical properties Based on the electronic ground state determined by HSE06, the complex frequency-dependent dielectric function e(x) = e1(x) + ie2(x) was investigated. The imaginary part e2(x) is determined by a summation over the transitions between the occupied and unoccupied states. The real part e1(x) is obtained by the Kramers–Kronig transformation [29,30]. Other optical spectrums such as refractive index n(x), absorption coefficient a(x) can be derived from e1(x) and e2(x) as follows [31]:

1 nðxÞ ¼ pffiffiffi 2

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 12 e21 ðxÞ þ e22 ðxÞ þ e1 ðxÞ

ð3Þ

Table 2 Calculated conduction band and valence band eigenvalues at high-symmetry points in the Brillouin zones of GaP and CdSiP2 using HSE06 and available experimental values. ZB

C1c X1c(z) X1c (xy) L1c

C15v a b

Ref. [24]. Ref. [25].

CH

C1c C3c T1c + T2c N1c

C4v C5v

GaP

CdSiP2

HSE06

Expt.a

HSE06

Expt.

2.56 2.41 2.41 2.43 0.00

2.89 2.39 2.39 2.64 0.00

2.79 2.06 2.70 3.29 0.00
0.27

0.00
0.20b

J. Xiao et al. / Computational Materials Science 117 (2016) 472–477

475

Fig. 4. PDOS of CdSiP2 calculated within HSE06. The Fermi level is set to zero.

 Fig. 6. PBE calculated ELF contour plot projected on the ð110Þ plane of a unit cell.

Fig. 5. Band decomposed charge density of (a) VBM and (b) CBM calculated within HSE06. The isosurfaces correspond to 0.025 e/Å3 for (a). The value of isosurfaces for (b) is set at 0.02 e/Å3 to visualize the orbital characters of CBM. The big red balls, media blue balls and small green balls represent Cd, Si and P atoms. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Table 3 HSE06 calculated Bader charge of CdSiP2. Atomic species

Cd

Si

P

Electronic charge (e)

11.41

2.20

6.19 6.20

12 pffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e21 ðxÞ þ e22 ðxÞ
e1 ðxÞ

aðxÞ ¼ 2x

ð4Þ

The dielectric functions along different directions of polarized light are presented in Fig. 7. Due to the symmetry, CdSiP2 has one optical axis being coincident with z-axis. So the dielectric tensor in the principal coordinate system has two independent components, ex = ey and ez. For ex = ey and ez, the electric fields E are perpendicular and parallel to the optical axis (c-axis), respectively. As the peaks originate from the contribution of transitions in the entire k-space, the transitions between the high-symmetry points

Fig. 7. (a) Real part and (b) imaginary part of the dielectric function for light polarized along different directions calculated using HSE06.

are not discussed here. There are two main peak complexes above the absorption edge in the spectrums of e2(x). The first peak complex between 4 eV and 4.6 eV almost derives from transitions between P-2p valence bands and the conduction bands. The peak complex B between 5.1 and 5.6 eV is mainly the contribution of

476

J. Xiao et al. / Computational Materials Science 117 (2016) 472–477

Fig. 8. Absorption coefficient for light polarized along different directions calculated using HSE06.

transitions from Si-3p and P-2p valence bands to the conduction bands. Another important thing is that our result shows obvious optical anisotropy in the region between 3.5 eV and 8 eV. This anisotropy derives from the transitions involve pz electrons or (px, py) electrons of P atoms and Si atoms. Now we turn to the static dielectric constants including both electronic and ionic contributions. The dielectric constants e1x ðx ! 0Þ and e1z ðx ! 0Þ are 8.460 and 8.246 separately. These two values represent the electronic contributions to the static dielectric constants. The ionic contributions to the static dielectric properties were obtained by density-functional perturbation theory [32]. Here the ionic contributions were calculated within GGA because GGA can provide a reliable result [33]. The ionic contributions values are 2.495 and 2.246 (E\c and Ekc). Thus, the total static dielectric constants are 10.955 and 10.492. Fig. 8 shows optical absorption coefficient. The shape and intensity differ along different polarization directions. The curves of optical absorption coefficient increase steadily at about 2.5 eV which is larger than the pseudodirect band gap calculated earlier. The possible reason for this phenomenon is that CdSiP2 has a pseudodirect band gap and transition from the VBM to CBM is symmetry forbidden as mentioned previously. Normally, the principal refractive indices of a uniaxial crystal are denoted nx = ny = no and nz = ne, where no and ne are referred as the ordinary and extraordinary refractive indices. The results of calculated refractive indices in Fig. 9 show that CdSiP2 crystal has a negative birefringence (Dn = ne
no < 0) below the absorption edge as expected. In practice, only the birefringence in the non-absorbing region (below the energy gap) is helpful in NLO application. However, there are discrepancies between calculated and experimental refractive indices. The primary source of discrepancies is the use of the random phase approximation (RPA) which neglects the electron–hole interaction to calculate the optical functions. The curves of the ordinary and extraordinary refractive indices (solid lines) shown in Fig. 9 were fitted by the Sellmeier equations. Normally, the Sellmeier equation for an IR NLO crystal is given by [34]

n2 ¼ A þ

B 1
C=k2

þ

D 1
E=k2

ð5Þ

Here, A represents the contribution from the electronic transitions above the band gap. B is oscillator strength and C is the square of the average absorption of the wavelength for the electronic transition. D and E represent similar as B and C, but for an ionic term [34,35]. The dispersion effects in the long-wavelength region are dominated by the ionic effect. Optical calculation using

Fig. 9. Refractive index versus wavelength for no and ne. The squares and triangles indicate no and ne computed by HSE06. The solid lines are fitted results by the Sellmeier equations. The dash and dot-dash lines are fitted by experimental dispersion equations from Ref. [2].

Table 4 Sellmeier coefficients derived from the calculated refractive indices. CdSiP2

A

B

C (lm2)

no ne

3.358 3.954

5.102 4.292

0.066 0.076

DFT above is based on the summation over the electronic transitions and doesn’t take into account the effects of lattice vibrations. Based on this, the calculated ordinary and extraordinary refractive indices were fitted by

n2 ¼ A þ

B 1
C=k2

ð6Þ

The Sellmeier coefficients A, B, C are presented in Table 4. The dispersion curves obtained by fitting calculated data are almost flat in the mid-IR and far-IR region, while the experimental curves have a declining trend in the same wavelength range. In order to match experimental dispersion curves, D is identified as the average phonon frequency, E is a fix value for all chalcopyrite compounds and A is corrected according to experimental indices of refraction in the method proposed by Lambrecht et al. [36]. But considering that this method is empirical and values of adjustment depend on experimental results, this adjustment is not included in our work.

4. Conclusions We have calculated the structural, electronic, bonding and optical properties within the HSE06 hybrid functional. Our results show that HSE06 with 25% short-range HF exchange can provide good predictions of equilibrium volume, bulk modulus and band gap. By comparing with the parent compound GaP, the correspondence between CBMs of GaP and CdSiP2 indicates that CdSiP2 has a pseudodirect band gap. The analyses of Bader charge and the electronic localization function indicate that CdSiP2 is not a pure ionic compound in which bonding types of P–Si and P–Cd are covalent and ionic, respectively. The observed features including main peaks and optical anisotropy in optical spectrums are mainly associated with p states. Through a systematic and accurate theoretical study, it will be helpful to clarify the puzzling issues of CdSiP2.

J. Xiao et al. / Computational Materials Science 117 (2016) 472–477

Acknowledgements This work was supported by the National Natural Science Foundation Programs of China (No. 51172149). Calculations were performed at the Center for High Performance Computing at Physics Discipline of SCU. References [1] V. Petrov, F. Noack, I. Tunchev, P. Schunemann, K. Zawilski, in: SPIE LASE: Lasers and Applications in Science and Engineering, International Society for Optics and Photonics, 2009, pp. 71970M-71970M-71978. [2] K.T. Zawilski, P.G. Schunemann, T.C. Pollak, D.E. Zelmon, N.C. Fernelius, F.K. Hopkins, J. Cryst. Growth 312 (2010) 1127–1132. [3] L. Fan, S. Zhu, B. Zhao, B. Chen, Z. He, H. Yang, G. Liu, X. Wang, J. Cryst. Growth 364 (2013) 62–66. [4] G. Zhang, H. Ruan, X. Zhang, S. Wang, X. Tao, CrystEngComm 15 (2013) 4255– 4260. [5] F. Chiker, B. Abbar, A. Tadjer, H. Aourag, B. Khelifa, Mater. Sci. Eng. 98 (2003) 81–88. [6] Z.-L. Lv, Y. Cheng, X.-R. Chen, G.-F. Ji, Comput. Mater. Sci. 77 (2013) 114–119. [7] Z. He, B. Zhao, S. Zhu, B. Chen, H. Hou, Y. Yu, L. Xie, Comput. Mater. Sci. 72 (2013) 26–31. [8] V. Shaposhnikov, A. Krivosheeva, V. Borisenko, J.-L. Lazzari, F.A. d’Avitaya, Phys. Rev. B 85 (2012) 205201. [9] Z. Wu, R. Cohen, D. Singh, Phys. Rev. B 70 (2004) 104112. [10] K. Hummer, J. Harl, G. Kresse, Phys. Rev. B 80 (2009) 115205. [11] G. Kresse, D. Joubert, Phys. Rev. B 59 (1999) 1758. [12] G. Kresse, J. Furthmüller, Comput. Mater. Sci. 6 (1996) 15–50.

477

[13] J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [14] J. Heyd, G.E. Scuseria, M. Ernzerhof, J. Chem. Phys. 118 (2003) 8207– 8215. [15] K. Momma, F. Izumi, J. Appl. Cryst. 44 (2011) 1272–1276. [16] A. Chahed, O. Benhelal, S. Laksari, B. Abbar, B. Bouhafs, N. Amrane, Physica B 367 (2005) 142–151. [17] F. Birch, Phys. Rev. 71 (1947) 809. [18] I. Hamdi, N. Meskini, Physica B 405 (2010) 2785–2794. [19] M. Böhm, G. Huber, A. MacKinnon, O. Madelung, A. Scharmann, E. Scharmer, New Serious, Group III 17 (1985). [20] M. Bettini, W. Holzapfel, Solid State Commun. 16 (1975) 27–30. [21] A. Kosobutsky, Y.M. Basalaev, J. Phys. Chem. Solids 71 (2010) 854–861. [22] S. Chen, X. Gong, A. Walsh, S.-H. Wei, Phys. Rev. B 79 (2009) 165211. [23] A. Shileika, Surf. Sci. 37 (1973) 730–747. [24] X. Zhu, S.G. Louie, Phys. Rev. B 43 (1991) 14142. [25] O. Madelung, U. Rössler, M. Schul, in, Springer, 2000. [26] G. Henkelman, A. Arnaldsson, H. Jónsson, Comput. Mater. Sci. 36 (2006) 354– 360. [27] A.D. Becke, K.E. Edgecombe, J. Chem. Phys. 92 (1990) 5397–5403. [28] L. De Santis, R. Resta, Surf. Sci. 450 (2000) 126–132. [29] M. Gajdoš, K. Hummer, G. Kresse, J. Furthmüller, F. Bechstedt, Phys. Rev. B 73 (2006) 045112. [30] C. Okoye, Physica B 337 (2003) 1–9. [31] V. Wang, W. Xiao, D.-M. Ma, R.-J. Liu, C.-M. Yang, J. Appl. Phys. 115 (2014) 043708. [32] S. Baroni, S. De Gironcoli, A. Dal Corso, P. Giannozzi, Phys. Rev. B 73 (2001) 515. [33] J. Paier, R. Asahi, A. Nagoya, G. Kresse, Phys. Rev. B 79 (2009) 115126. [34] G.C. Bhar, J. Phys. D Appl. Phys. 13 (1980) 455. [35] R. Xiao, H. Wu, Y. Ni, Z. Wang, C. Huang, M. Qi, C. Ge, J. Appl. Phys. 117 (2015) 135702. [36] W.R. Lambrecht, X. Jiang, Phys. Rev. B 70 (2004) 045204.