Computational Materials Science 106 (2015) 155–160
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
YN nanostructure formation on the GaN(0 0 0 1) surface: First principles studies J. Guerrero-Sánchez a,⇑, Gregorio H. Cocoletzi a, J.F. Rivas-Silva a, Noboru Takeuchi b a b
Benemérita Universidad Autónoma de Puebla, Instituto de Física ‘‘Ing Luis Rivera Terrazas’’, Apartado Postal J-48, Puebla 72570, Mexico Universidad Nacional Autónoma de México, Centro de Nanociencias y Nanotecnología, Apartado Postal 2681, Ensenada, Baja California 22800, Mexico
a r t i c l e
i n f o
Article history: Received 2 December 2014 Received in revised form 29 April 2015 Accepted 30 April 2015 Available online 22 May 2015 Keywords: Surfaces Yttrium nitride Gallium nitride Adsorption Diffusion Surface formation energy
a b s t r a c t We have investigated the Y adsorption and YN thin film formation on the GaN(0 0 0 1)-2 2 surfaces using first principles total energy calculations within the density functional theory. Ga-rich conditions are modeled with a Ga-bilayer terminated GaN(0 0 0 1)-2 2 surface. N-rich, and intermediate growth conditions were studied using an ideally GaN bilayer terminated GaN(0 0 0 1)-2 2 surface. On the ideally terminated surface, when the Y atom is constrained on top of the surface, the Y on a second-layer nitrogen atom (T4-site) is the most favorable structure. However, when the Y atom migrates into the first Ga monolayer, the displaced Ga ad-atom occupies a site bonded with three Ga atoms of the first monolayer (T4-2-site) as the most stable structure. Under Ga-rich conditions the Y atom occupies a site on top of a second-monolayer Ga (T4-site) as the most stable structure. Nevertheless, in the energetically most favorable configuration, the Y atom replaces a third inner layer Ga atom bonding to N atoms and forming an YN pair. When the coverage is increased to a full Y monolayer, three different YN configurations are formed. Under Ga-rich conditions an YN bilayer is formed underneath a Ga bilayer, in the intermediate conditions an YN bilayer is formed under a Ga-T4 layer, and under N-rich conditions a w-YN bilayer is formed on top of the ideal surface. The stable configurations density of states shows that the metallic characteristic is preserved. Ó 2015 Elsevier B.V. All rights reserved.
1. Introduction Wurtzite phase gallium nitride (GaN) is a semiconductor with a wide band gap of 3.4 eV, large bulk modulus and high thermal conductivity; it is an excellent material for technological applications. GaN is currently found in a large variety of applications such as high power, high temperature and high electronic mobility devices [1–6]. Also, a strong inter-sub-band absorption peak located at 675 meV in AlGaN/GaN super-lattices has been found, so that it is possible to apply these heterostructures in devices which work in the mid and near-infrared [7,8]. Moreover, gallium nitride may form alloys with other group IIIA elements, especially with indium, resulting in a material with continuous variation of energy gap from 0.7 eV to 3.4 eV, this band gap engineering can be applied in high-efficiency photovoltaic devices [9]. On the other hand, transition metals have been proposed to substitute group IIIA elements on these alloys with GaN. In particular gap engineering has been realized in ScxGa1xN alloys where the band gap can be varied linearly in the range of 2–3.4 eV with variations in the Sc ⇑ Corresponding author. E-mail address:
[email protected] (J. Guerrero-Sánchez). http://dx.doi.org/10.1016/j.commatsci.2015.04.050 0927-0256/Ó 2015 Elsevier B.V. All rights reserved.
concentration [10]. First principles calculations have shown that ScN has a rock-salt structure [11,12]. But, it forms a hexagonal structure in ScGaN alloys [13,14]. Additionally, it is possible to reduce defects in heteroepitaxial non-polar and semi-polar GaN films with wurtzite-like interlayers of ScN [15]. A recent report has appeared [16] that deals with the ScN nanostructure formation on the GaN(0 0 0 1) surface under both Ga-rich and N-rich conditions. The first principles total energy calculations showed that Sc atoms prefer to be in the GaN sub-surface, exchanging places with the first gallium monolayer to form ScN bilayers [16]. It is known that YN is a transition metal compound with semiconductor behavior similar to ScN [17]. Several groups have investigated YN using first principles calculations with results indicating that YN crystallizes in a rock-salt structure with a metastable wurtzite structure [17–19]. Recently YN thin films have been fabricated by reactive laser ablation [20]. The lattice parameter mismatch between YN and GaN is of the order of 8%, similar to that of InN and GaN. Since InGaN alloys have been already grown and proposed for applications in photovoltaic devices [9], it is expected that YN may have a similar behavior when alloying or forming heterostructures with GaN. Moreover, the increase interest in those YGaN and YN/GaN systems has motivated first principles
156
J. Guerrero-Sánchez et al. / Computational Materials Science 106 (2015) 155–160
studies, with the studies showing that GaxY1xN alloys display a stable wurzite-like structure between moderate and high Ga concentrations [21]. It has been also reported a modulation in the band gap only by increasing the Ga content [21]. It has been showed that 4d transition-metal doped GaN nanotubes exhibit an important structural distortion at the impurity vicinities, in particular; Y-doped GaN nanotubes have a semiconductor behavior [22]. Taking into account these interesting findings we study the Y adsorption and the initial stages of the YN structure formation on the GaN(0 0 0 1) surface using first principles total energy calculations. The paper is organized as follows: In Section 2 we describe the computational procedure, in Section 3 we present and discuss results, and in Section 4 we make conclusions.
2. Computational details In this work we perform calculations using the density functional theory [23,24], as implemented in the PWscf code of the Quantum ESPRESSO package [25]. We have used Vanderbilt ultra-soft pseudopotentials to deal with the electron–ion potential [26], and the exchange–correlation interactions are treated according to the generalized gradient approximation as implemented in a Perdew–Burke–Ernzerhof (PBE) functional [27]. The electron Kohn–Sham states are expanded in plane waves with a kinetic-energy cutoff equal to 30 Ry and for the charge density 240 Ry. An equally-spaced k-point grid, as proposed by Monkhorst–Pack of 3 3 1 is used to sample the irreducible Brillouin zone [28]. Energy and force convergence criteria have been applied to the structural optimizations; the total energy changes are less than 1 mRy and the Hellmann–Feynman force is less than 0.002 Ry/Å. On the other hand, a smearing [29] width of 0.04 Ry has been used. All these parameters allow us to have a good confidence to compare energy differences in the models. Our studies consider the previously calculated and reported parameters for GaN in the wurtzite phase [16]; a = 3.22 Å, c/a = 1.63 and u = 0.38, which are in good agreement with the experimental data [30]. In order to obtain the density of states (DOS) and the projected DOS we increase the k-points grid to 21 21 1, also a Methfessel–Paxton smearing [29] with a smearing width of 0.07 eV was used for the Brillouin-zone integrations. Studies are done using a GaN(0 0 0 1)-2 2 surface model. The supercell slab consists of four GaN bilayers, with four atoms per layer. Under Ga-rich conditions a bilayer of Ga is added on top of the ideally terminated surface. We test slabs with 4, 5, 6, and 7 bilayers, after-relaxation results show that in the third, fourth, fifth and sixth bilayer, respectively, the average change in Ga–N distances is 0.001 Å, then it is unnecessary to use slabs with large number of bilayers for calculations. The lowest N-layer of the slab is passivated with pseudo-hydrogen atoms with a fractional charge of 0.75e. To determine the suitable vacuum space, we have tested different vacuum gaps 11, 13, and 15 Å. Results show that the energy differences between models (Y on GaN) are the same as those reported in the section of results. Then, it is only required an empty space of 11 Å to avoid the undesirable charge transfer between adjacent slabs. Also, the vacuum gap is enough to neglect the dipole correction, provided we have almost the same empty gap as used in [31]. In order to simulate a bulk-like behavior we freeze the bottom GaN bilayer and the pseudo-hydrogen atoms. The ad-atoms binding energies, at stable and metastable sites, have been calculated after fully relaxing the atomic positions. For unstable sites, we keep fixed the lateral positions of the ad-atoms at a given position, and allowing all other coordinates to relax. In this way we map the potential energy surface (PES). The minimum values of the PES give the adsorption sites, while maximum values and saddle sites give the diffusion barriers.
To determine the energetically most stable configurations with different number of atoms, we have computed the surface formation energies as reported in [32]. Our calculations are performed in thermal equilibrium, considering that lGaN(bulk) = lGa + lN, with lGa 6 lGaðbulkÞ and lN 6 lN2 ðmoleculeÞ in order to avoid Ga and N precipitations. In this way we may write,
lN2 ðmoleculeÞ
DHGaN , f
where
DHGaN f
lGaNðbulkÞ ¼ lGaðbulkÞ þ
= 0.99 eV is the calculated GaN
heat of formation. At T = 0 we approximate G(lY, lGa, lN) E(lY, lGa, lN) to write the relative formation energy as:
Eformation ¼ Etot Eref DnY lY DnGa lGa DnN lN where Etot is the total energy of the configuration under study, Eref is the total energy of the reference system, in our case is the clean GaN(0 0 0 1) surface. The surface formation energies for all configurations are calculated with lY = lY(bulk) and the limits
lGaðbulkÞ DHGaN 6 lGa 6 lGaðbulkÞ which correspond to the interval f between N-rich conditions (lower limit) and Ga-rich conditions (upper limit). 3. Results In this section we present results of the adsorption of a single Y ad-atom and its diffusion when is constrained on top of the GaN(0 0 0 1) surface, the possibility of migration into inner sites, the effects of increasing the Y coverage, and the initial stages of the YN bilayer growth under Ga-rich, moderately Ga-rich and N-rich conditions. 3.1. Y adsorption on the GaN bilayer terminated GaN(0 0 0 1) surface It has been demonstrated that under N-rich growth conditions the GaN(0 0 0 1) surface exhibits a 2 2-N ad-atom structure, while in moderately Ga-rich conditions a 2 2-Ga ad-atom model is the more favorable [33]. Therefore, the GaN bilayer terminated surface represents a very realistic model of the surface during growth under N-rich and moderately Ga-rich conditions. Fig. 1A shows a schematic representation of the surface model with indications of the high symmetry adsorption sites and the diffusion pathway in dotted lines. We first study the adsorption of Y ad-atoms at high symmetry sites: T4, H3, Top, and Bridge (Br). It is found that after relaxation, Y adsorption on the T4 site (forming bonds with three first layer Ga atoms) is the most favorable, as it has been found with other ad-atoms [34]. The relaxed H3 site is metastable, while the Top is an unstable site. When compared with the total energy of the Y ad-atom on a T4 site (from now on, this energy is used as reference) these sites have higher energies of 0.15 eV and 1.61 eV, respectively (see Table 1). To study the Y ad-atom surface diffusion, we have calculated the PES profile along the dotted line in Fig. 1A. An Y ad-atom diffusing from T4 to H3 encounters an energy barrier of 0.51 eV at the Br site as shown in Fig. 2A, the path between T4 and H3 was divided in 5 steps to precise the position of the Bridge site. Diffusion along a path through the Top site is highly non probable because of the large potential barrier of 1.61 eV of height. Now, we consider the Y ad-atom migration to the first Ga monolayer by replacing a Ga atom. The calculated YN heat of formation (3.00 eV) is larger than the heat of formation of GaN (0.99 eV), making the formation of YN instead of GaN thermodynamically more favorable, and therefore giving rise to the Y ad-atom migration to inner GaN layers. The replacing of a first layer Ga atom by an Y atom results in new different T4, H3, Top and Bridge sites. We show them in Fig. 1B and label them in two groups: T4-1, H3-1, Top-1, and Br-1 and T4-2, H3-2, Top-2, Br-2,
157
J. Guerrero-Sánchez et al. / Computational Materials Science 106 (2015) 155–160
Fig. 1. GaN(0 0 0 1) surface (A) GaN bilayer-terminated surface, and (B) Ga is replaced by an Y atom on the GaN bilayer-terminated surface. (C) Ga bilayer terminated surface. The dotted lines in all cases (top view) indicate the possible ad-atom diffusion pathway.
Table 1 Relative energies (in eV) of Y ad-atom at high symmetry sites on a GaN bilayer terminated GaN(0 0 0 1) surface, the reference energy corresponds to the Y-T4 site over the surface. Site
Y above
Site
Y below, Ga ad-atom bonded with the Y atom
Site
Y below, Ga ad-atom bonded only with Ga
T4 H3 Top Bridge
0.00 0.15 1.61 0.51
T4-1 H3-1 Top-1 Bridge-1
2.21 2.14 0.89 1.80
T4-2 H3-2 Top-2 Bridge-2, Bridge-3
3.44 2.89 1.84 2.62, 2.95
and Br-3. In the first group, the displaced Ga ad-atom bonds to the Y and first layer Ga atoms, while in the second group it bonds with first layer Ga atoms only. When the Y ad-atom exchanges sites with a Ga atom, it is energetically more favorable for the displaced
Ga atom to occupy the second kind of high symmetry sites. Moreover each of these configurations has a lower energy than that of the Y at the T4 site. Results show that the most stable site is the T4-2, as summarized in Table 1.
158
J. Guerrero-Sánchez et al. / Computational Materials Science 106 (2015) 155–160
Fig. 2. Diffusion pathways (in eV) for an Y ad-atom at high symmetry sites on (A) a GaN bilayer terminated GaN(0 0 0 1) surface and (B) a double Ga layer on GaN bilayer terminated GaN(0 0 0 1) surface. The dots of the high symmetry sites were plotted twice in order to show the diffusion path and the line is only a guide of the path followed by the ad-atom at the high symmetry sites (H3-Br-T4-Br-H3).
3.2. Y adsorption on the Ga bilayer terminated GaN(0 0 0 1) surface Reports have shown that in the GaN growth under Ga-rich conditions the result is a smooth surface which is due to the formation of a fluid-like Ga double monolayer at the growth front that acts as surfactant, lowering the surface energy. It has been shown that the GaN(0 0 0 1) surface displays a pseudo-1 1 periodicity and surface corrugation as indicated by the scanning tunneling microscopy (STM) images [35]. Models have been proposed to explain the experimental findings: a fluidic contracted Ga bilayer at the surface [35], the topmost layer to be biaxially contracted, while the layer beneath remains bulk-like [36], and a uniaxially contracted configuration [37]. Therefore, taking into account these findings we study the adsorption under Ga-rich conditions, we have considered the adsorption and diffusion of Y ad-atoms on a Ga-bilayer terminated GaN(0 0 0 1) surface using a 2 2 unit cell as shown in Fig. 1C. The relaxed structures display similar atomic configurations to those obtained on the ideally terminated surface: the Y ad-atoms are adsorbed at T4 sites in the most stable geometry. The H3 and Bridge sites are metastable, with total energies of 0.06 eV and
0.07 eV higher than the T4 geometry, respectively. Since the energy barrier is small, the Y ad-atom diffusion is highly favorable, see Fig. 2B. To investigate the Y migration into the inner layers, we substitute a Ga atom at inner layers by an Y atom, with the displaced Ga atom on top of the surface. Results show that the most stable configuration corresponds to the Y atom replacing a third layer atom and the displaced Ga atom on a T4 site. In this configuration the Y atom forms YN bonds with three N atoms. Y at first and second layers result in higher energies as indicated in Table 2.
Table 2 Relative energies (in eV) of Y ad-atom at high symmetry sites on a Ga bilayer terminated GaN(0 0 0 1) surface, the reference energy corresponds to the Y-T4 site over this surface. Site
Y above
Layer
Site
Y at inner layers
T4 H3 Bridge
0.00 0.06 0.07
First Second Third
T4(2) H3 T4
1.80 1.10 2.16
J. Guerrero-Sánchez et al. / Computational Materials Science 106 (2015) 155–160
3.3. Surface formation energies To investigate the Y coverage effects under nitrogen rich conditions and moderately gallium rich conditions, we consider the adsorption of ½, 3=4 and 1 monolayer (ML) of Y on a GaN bilayer terminated surface. It is found that 1 ML of Y deposited at the T4 site
159
yields the most favorable structure, with total energies of 0.24 eV/(1 1) and 0.28 eV/(1 1) lower than these of 3=4 ML and ½ ML coverages, respectively. However, it is found that a configuration with all Y atoms replacing the first Ga layer and the displaced atoms occupying T4 sites is energetically more favorable. This configuration results in the formation of an YN bilayer under a Ga monolayer, and as seen in Fig. 3, it is stable for moderately Ga-rich conditions, in the interval between 0.03 and 0.4 eV of Ga chemical potential. Under nitrogen rich conditions, there is another possibility for YN growth, in which a N monolayer occupies top sites above Ga atoms, followed by an Y layer, resulting in the formation of an YN bilayer on top (this is equivalent to the replacement of the topmost Ga layer). As seen in Fig. 3, this configuration is the most favorable structure under nitrogen rich conditions. We have also investigated the growth of Y and YN thin films under Ga rich conditions, by considering the adsorption of ½, 3=4 and 1 monolayer (ML) of Y on a Ga bilayer terminated surface. Under these conditions, the formation of a full monolayer of Y is also more stable than adsorption of ½ and 3=4 ML of Y by 0.06 eV/(1 1) and 0.09 eV/(1 1), respectively. However it is found that a configuration in which the third Ga layer is replaced by the Y layer is the most stable configuration under gallium rich conditions as shown in Fig. 3. In this structure, there is a formation of an YN bilayer underneath a Ga bilayer with the topmost Ga atoms occupying H3 sites. 3.4. Electronic properties
Fig. 3. Surface formation energies calculated as a function of the Ga chemical potential for the GaN(0 0 0 1) surface. The energy for the ideal GaN bilayer terminated surface is used as zero.
In this section we present the total density of states (DOS) and partial density of states (PDOS) of the most favorable atomic configurations to describe the electronic structure changes induced by the YN bilayer formation on GaN. In Fig. 4 we plot the DOS and PDOS of the YN bilayer formation: (A) on the GaN bilayer terminated surface, (B) under a Ga-T4 layer, and (C) under a Ga bilayer. In (A) a metallic behavior is obtained, which is mostly induced by the Y atoms dangling bonds on top of the structure. The main contribution to the occupied states comes from the N-p orbitals, whereas the main contributions to the unoccupied band are due to Y-d and Ga-p orbitals. The total DOS of geometry (B), shown in Fig. 4B, exhibits also metallic characteristics mostly produced by the Ga-p orbitals at the vicinity of the Fermi level. Similar
Fig. 4. Density of states (DOS) and Partial density of states (PDOS) of the (A) YN bilayer on GaN(0 0 0 1) surface, (B) YN under a Ga monolayer, and (C) YN under a Ga bilayer (Ga-rich conditions). The Fermi level is set at zero eV.
160
J. Guerrero-Sánchez et al. / Computational Materials Science 106 (2015) 155–160
to configuration (A) the principal contribution to the occupied states comes from the N-p orbitals, and for the unoccupied states the Ga-p and Y-d orbitals are those with largest contributions. Finally in Fig. 4C, we show the density of states of configuration (C). This system exhibits also metallic characteristics; a large population of Ga-p states exists near the Fermi level, as induced by the extra Ga atoms. At the unoccupied orbitals the Ga-p and Y-d states make the most important contributions. On the other hand, N-p orbitals also make large contribution to the occupied states. In order to explain the stability of the w-YN bilayer configuration (Fig. 4A) we have analyzed the partial density of states, it is clear that there exist a hybridization of the Y-d and N-p orbitals in the energy range from 2.5 to 5 eV. This fact supports the Y–N bonding formation, a similar behavior is obtained for (B) and (C) configurations. As expected, the YN formation on GaN induce changes in the DOS, showing that in all DOS the Y-d and Ga-p orbitals make the most important contributions to the empty states, whereas N-p orbitals make the largest contribution to the occupied states. 4. Conclusions We have presented first principles total energy calculations to investigate the Y and YN film formation on the GaN(0 0 0 1) surface. We have studied the single Y atom adsorption on an ideally GaN bilayer terminated GaN(0 0 0 1) surface considering N rich, and on a Ga bilayer terminated GaN(0 0 0 1) surface, corresponding to Ga rich conditions. In both cases the most stable configurations correspond to the Y adsorption at the T4 site. However the Y migration to inner layers results energetically more favorable in both cases when the Y bonds to N atoms to form YN. A completion of a full monolayer is always more stable than the adsorption of ½ or 3=4 monolayers. Under these conditions three different configurations are possible: a w-YN bilayer terminated configuration for N rich conditions, an YN bilayer under a Ga-T4 layer for the intermediate case and an YN bilayer underneath a Ga bilayer for Ga rich conditions. In all three geometries, the YN bilayers are in the wurtzite-like phase with the surfaces showing metallic characteristics. Acknowledgements G.H.C. acknowledges the financial support of VIEP-BUAP, grant 31/EXC/06-G, CONACYT project #223180 and Cuerpo Académico Física Computacional de la Materia Condensada (BUAP-CA-191). NT thanks DGAPA project IN103512-3 and Conacyt Project 164485 for partial financial support. JGS thanks M.C. Francisco Sánchez-Ochoa for useful discussions. Calculations were performed in the DGCTIC-UNAM supercomputing center, the computer center of the Instituto de Física BUAP, and the Centro Nacional de Supercomputo, IPICYT.
References [1] Wataru Saito, Yoshiharu Takada, Masahiko Kuraguchi, Kunio Tsuda, Ichiro Omura, Tsuneo Ogura, Hiromichi Ohashi, IEEE Trans. Electron Devices 50 (2003) 2528. [2] Liang Pang, Kyekyoon Kim, J. Phys. D Appl. Phys. 45 (2012) 045105. [3] Hong-An Shih, M. Kudo, M. Akabori, Toshi-kazu Suzuki, Jpn. J. Appl. Phys. 51 (2012) 02BF01. [4] Yuji Wang, Wu Lu, Phys. Status Solidi A 208 (2011) 1623–1625. [5] Chi-Shiang Hsu, Huey-Ing Chen, Chung-Fu Chang, Tai-You Chen, Chien-Chang Huang, Po-Cheng Chou, Wen-Chau Liu, Sens. Actuat. B 165 (2012) 19–23. [6] A. Zanandrea, E. Bahat-Treidel, F. Rampazzo, A. Stocco, M. Meneghini, E. Zanoni, O. Hilt, P. Ivo, J. Wuerfl, G. Meneghesso, Microelectron. Reliab. 52 (2012) 2426. [7] C. Bayram, J. Appl. Phys. 111 (2012) 013514. [8] C. Edmunds, L. Tang, D. Li, M. Cervantes, G. Gardner, T. Paskova, M.J. Manfra, O. Malis, J. Electron. Mater. 41 (2012) 881. [9] Dirk V.P. McLaughlin, Joshua M. Pearce, Metall. Mater. Trans. A 44A (2013) 1947. [10] M.E. Little, M.E. Kordesch, Appl. Phys. Lett. 78 (2001) 2891. [11] Noboru Takeuchi, Phys. Rev. B 65 (2002) 045204. [12] R. Mohammad, S. Katircioglu, Condens. Matter Phys. 14 (2011) 23701. [13] Costel Constantin, Muhammad B. Haider, David Ingram, Arthur R. Smith, Nancy Sandler, Kai Sun, Pablo Ordejón, J. Appl. Phys. 98 (2005) 123501. [14] Maria Guadalupe Moreno-Armenta, Luis Mancera, Noboru Takeuchi, Phys. Status Solidi (b) 238 (2003) 127–135. [15] M.A. Moram, C.F. Johnston, M.J. Kappers, C.J. Humphreys, J. Cryst. Growth 311 (2009) 3239–3242. [16] J. Guerrero-Sánchez, Gregorio H. Cocoletzi, J.F. Rivas-Silva, Noboru Takeuchi, Appl. Surf. Sci. 268 (2013) 16–21. [17] S. Zerroug, F. Ali Sahraoui, N. Bouarissa, Appl. Phys. A 97 (2009) 345–350. [18] G. Soto, M.G. Moreno-Armenta, A. Reyes-Serrato, Comput. Mater. Sci. 42 (2008) 8–13. [19] Luis Mancera, Jairo A. Rodríguez, Noboru Takeuchi, J. Phys.: Condens. Matter 15 (2003) 2625–2633. [20] W. De La Cruz, J.A. Díaz, L. Mancera, N. Takeuchi, G. Soto, J. Phys. Chem. Solids 64 (2003) 2273–2279. [21] Youcef Cherchab, Mohamed Azzouz, Rafael González-Herná¡ndez, Khedija Talbi, Comput. Mater. Sci. 95 (2014) 509–516. [22] Seyyed Hassan Saberi, Seyyed Mahdy Baizaee, Hamideh Kahnouji, Superlattices Microstruct. 74 (2014) 52–60. [23] U. von Barth, Phys. Scripta T109 (2004) 9–39. [24] M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias, J.D. Joannopoulos, Rev. Mod. Phys. 64 (1992) 1045. [25] Paolo Giannozzi et al., J. Phys.: Condens. Matter 21 (2009) 395502. [26] David. Vanderbilt, Phys. Rev. B 11 (1990) 7892. [27] J.P. Perdew, S. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. [28] H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. [29] M. Methfessel, A.T. Paxton, Phys. Rev. B 40 (1989) 3616. [30] H.P. Maruska, J.J. Tietjen, Appl. Phys. Lett. 15 (1969) 327. [31] Lennart Bengtsson, Phys. Rev. B 59 (1999) 301. [32] Guo-Xin Qian, R.M. Martin, D.J. Chadi, Phys. Rev. B 38 (1988) 7649. [33] A.R. Smith, R.M. Feenstra, D.W. Greve, J. Neugebauer, J.E. Northrup, Phys. Rev. Lett. 79 (1997) 3934. [34] N. Takeuchi, A. Selloni, T.H. Myers, A. Doolittle, Phys. Rev. B 72 (2005) 115307. [35] A.R. Smith, R.M. Feenstra, D.W. Greve, M.S. Shin, M. Skowronski, J. Neugebauer, J.E. Northrup, J. Vac. Sci. Technol. B 16 (1998) 2242. [36] J.E. Northrup, J. Neugebauer, R.M. Feenstra, A.R. Smith, Phys. Rev. B 61 (2000) 9932. [37] J.A. Rinehimer, M. Widom, J.E. Northrup, R.M. Feenstra, Phys. Status Solidi (b) 245 (2008) 920.