- Email: [email protected]
Comparative study of the empirical interatomic potentials and density-functional simulations of divacancy and hexavacancy in silicon
Physica B 406 (2011) 467–470
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Comparative study of the empirical interatomic potentials and density-functional simulations of divacancy and hexavacancy in silicon Chaoying Wang a,n, Zhenqing Wang a, Qingyuan Meng b a b
College of Aerospace and Civil Engineering, Harbin Engineering University, Harbin 150001, PR China Department of Astronautical Science & Mechanics, Harbin Institute of Technology, Harbin 150001, PR China
a r t i c l e in f o
abstract
Article history: Received 8 April 2010 Received in revised form 4 October 2010 Accepted 5 November 2010
The comparison between three classical potentials and density functional theory (DFT) is performed mainly for divacancy and hexavacancy in Si crystal. According to their performances on the formation energies and structural properties (distortion magnitude, relaxation volume and symmetry), the limitations and validities of classical potentials are discussed. It is found that the outward relaxation directions of EDIP and Tersoff (T3) are contrary to the DFT and Stillinger–Weber (SW) directions (inward), which restrict their application in the structural property calculations. Except for the divacancy symmetries, the results of SW are in agreement with the DFT results. Thus, it can be concluded that SW should be the best potential to describe V2 and V6. & 2010 Elsevier B.V. All rights reserved.
Keywords: Vacancies DFT Empirical potential
1. Introduction Vacancies, such as divacancy (V2) and ring hexavacancy (V6), are fundamentally and technologically important defects in silicon. V2 is particularly attractive from the experimental point of view because it can be easily produced by electron irradiation, and moreover, it is quite stable. Nowadays, the structural properties or formation energies of V2 with a variety of charge states have been explored in the experiments [1,2] and ab initio calculations [3–6]. V6 is the smallest ring vacancy cluster. It is likely to be an efficient gettering center for a range of impurities, such as carbon, oxygen and copper atoms [7]. In addition, it can also be the precursors of extended defects, such as prismatic dislocations [8], H-related platelets [9] or photoluminescent (PL) center [10]. In the pioneer work, Chadi and Chang [11] used the dangling bonds counting (DBC) model and predicted that V6 should be more stable than other types of vacancy clusters. This was verified in the subsequent density functional theory (DFT) [12–14] and density functional based tight binding (DFTB) [15] calculations. This type of simple model of V6 is important since it can give a rough idea about how energetically favorable structures may look like. Moreover, the most important vacancy aggregates are those that are the most stable, as they are the most likely ones to survive high-temperature treatments. On the other hand, it is noted that the theoretical interest of vacancies is mainly focused on the DFT calculations. However, this method is unable to deal with the large-scale vacancies-related
problems, such as the interactions of vacancies with extended defects (grain boundary, dislocation [16,17] and stacking fault), self-interstitial migration in large systems and so on, due to large atomic system and long time period calculations. A possible alternative for these cases might be the use of empirical potentials, which are computationally much less expensive. However, the accuracy of those potentials cannot be ensured. Thus, for a particular defect, it is necessary to test the validity before their applications. For this purpose, a lot of tests has been done previously, e.g., dislocations [18,19], elastic constants [18,20], surfaces [20,21], thermo-mechanical properties [22] and so on. Moreover, Godet et al. [23] compared the Stillinger–Weber (SW) [24], environment-dependent interatomic potential (EDIP) [25] and Tersoff (T3) [26] with ab initio method under large shear. For vacancies, monovacancy (V1) is widely used for testing the ability of empirical potentials [18,20]. But for V2 and V6, to our knowledge, little work has been done on the determination of the best interatomic potentials. In this work, the structural properties and formation energies of V2 and V6 are calculated by the use of DFT and other three classical potentials. The limitations and validities of empirical potentials for modeling vacancies are discussed by means of the comparison of the calculated results. For the comparison between the DFT and empirical potentials results, only the neutral vacancies are considered in the present work.
2. Computational methods n
Corresponding author. Tel.: + 86 045182519956. E-mail address: [email protected] (C. Wang).
0921-4526/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2010.11.014
The first-principle calculations of the vacancies are conducted using SIESTA code [27,28]. The Perdew, Burke and Ernzerhof (PBE)
468
C. Wang et al. / Physica B 406 (2011) 467–470
[29] formulation of the generalized gradient approximation (GGA) is used for the electronic exchange and correlation. Core electrons are replaced by the non-local norm-conserving Troullier–Martins pseudopotential [30], which is factorized in the Kleinman– Bylander form [31]. Due to the requirement of more precise results, a real-space grid equivalent to 160 Ry plane wave cut-off is used. For the Brillouin zone (BZ) sampling, Monkhorst–Pack (MP) [32] mesh q3 shifted by (k0, k0, k0) is utilized. In all cases, atomic geometries are relaxed using the conjugate-gradient (CG) [33] ˚ method until interatomic forces are less than 0.01 eV/A. During the DFT calculations, the single-zeta (SZ), double-zeta (DZ) and doubled-zeta with polarization (DZP) basis sets can be used for Si atoms. The lattice constants (a0) with different basis sets ˚ respectively, which are are converged to 5.643, 5.556 and 5.494 A, 3.8%, 2.3% and 1.2%, within the experimental value 5.43 A˚ [34]. To remove spurious electronic interactions, the periodic supercell method is employed in the DFT calculations. With the different supercell sizes and a0 obtained, the total energies per atom with respect to q are plotted in Fig. 1. For all cases, q varies from 1 to 4, which corresponds to k0 ¼0 for q¼1 and k0 ¼0.5 for q¼2, 3, 4. Fig. 1 shows that the energy converges to different values with corresponding basis set. This implies that, except for computational techniques adopted, different a0 may lead to diverse results. In this work, the energy converges to 0.9 meV (DZ) and 0.04 meV (DZP) at q¼2 for 216 supercell. In addition, the mesh of 23 corresponds to a density of 0.031 A˚ 1, which is sufficient to make the total energy converge eventually [35]. In view of previous successful DFT calculations of V2 [4] and V6 [12,13], the 216 atoms supercell can be well converged eventually. Thus, the DFT calculations using the 216 supercell and DZP basis set at q¼2 are conducted in this work. Because of the great technological importance of Si, there are more than 30 empirical potentials, which have been developed and applied to different situations. In consideration to the advantages and limitations of each potential, the widely used and tested SW and T3 potentials [24] were applied for this work. Moreover, EDIP is known to describe both 301- and 901-partial dislocations in the glide set {1 1 1} [18], which was also used in this work. In fact, the SW potential is fitted to experimental properties of the diamondcubic (dc) and molten phases of silicon. Its formulation consists of a linear combination of two- and three-body terms and contains eight parameters. The EDIP also includes two- and three-body terms. which depend on the local atomic environment given by an effective coordination number. Its thirteen parameters in the functional form are obtained by fitting to a set of ab initio results,
Fig. 1. Convergence tests for the energy per atom as a function of Brillouin zone sampling for different sized supercells and basis sets.
which include various bulk phases and defect structures. The Tersoff potential consists of many-body interactions included in a bond order term. In addition, the eleven parameters are fitted to ab initio results for several Si polytypes. It has three versions in the history and the last version (T3) is used in this work. The total energies are calculated by the molecular dynamics (MD) method with the NPT (the atoms number, pressure and temperature are kept constant) ensemble and a time step of 1 fs. The vacancy is formed in the middle of MD model, which is 9.9 nm along the three vectors [1 0 0], [0 1 0] and [0 0 1]. In the calculations, the periodic boundary conditions (PBC) are applied in the three directions of the model. The total energy is minimized during the annealing process. First, atomic structure is relaxed for 10,000 steps at 0 K. And then, the temperature becomes 10 K higher at every 5000 steps until a desired temperature is achieved. By the reverse process, the temperature drops to 0 K and the energy is obtained. In the present work, zero-temperature formation energies of vacancies are computed using the expression Efn ¼ EðNnÞ
Nn EðNÞ N
ð1Þ
where E(N) is the total energy of the perfect crystal containing N atoms, E(N n) is the total energy of the crystal containing n vacant sites. As a measure of the variation in vacancy volumes resulting from the relaxation, we also calculate the relative relaxation volume
DO ¼ ðOO0 Þ=O0 100%
ð2Þ
where O0 and O are the volumes of the vacancy before and after relaxation, respectively.
3. Results and discussion 3.1. Divacancy Removing two adjacent atoms in the perfect crystal generates an unrelaxed V2, where atomic structures are of D3d symmetry. From Fig. 2(a), we can see that there are two sets of three distances between the three nearest-neighbors of each vacant site. During the DFT relaxation process, the Jahn–Teller (JT) distortion reduces the ideal D3d symmetry to C2h. Moreover, the pairing of the atoms 2–3 and 4–5 forms the large pairing (LP) [1] structure, where d2,3 ¼ d4,5 od1,2 ¼d1,3 ¼d4,6 ¼d5,6. During this process, the volume of V2 reduces by 38%. Except for the LP structure, Iwata et al. [3] claimed that the resonant bond (RB) structure was more stable than the DFT calculations in the supercell ranging from 64 to 1000 atoms. The RB structure is signaled by d2,3 ¼ d4,5 4d1,2 ¼d1,3 ¼d4,6 ¼d5,6. But in another large scale (1000-atom) DFT calculation [4], Wixom and Wright found that the LP structure was more stable than the RB. This conclusion was proved by other DFT calculations [5,6] and experiment [2]. The six distortion amplitudes relative to the unrelaxed
Fig. 2. Schematic view of (a) divacancy (V2), (b) hexavacancy (V6) and their nearestneighbor atoms.
C. Wang et al. / Physica B 406 (2011) 467–470
Table 1 Structural properties and formation energies (Ef2) of relaxed divacancy. Dd (%) indicates the fractional distortion magnitudes, which is relative to the ideal distance ˚ The signs and + indicate the inward and outward relaxation, respectively. 3.885 A.
469
Table 2 Structural properties and formation energies (Ef6) of relaxed hexavacancy. Dd (%) indicates the fractional distortion magnitudes, which is relative to the ideal distance ˚ The signs and + indicate the inward and outward relaxation, respectively. 3.885 A.
Method
Ef2 (eV)
Dd (%)
DO (%)
Symmetry
Method
Ef6 (eV)
Dd (%)
DO (%)
Symmetry
DFT SW EDIP T3
5.59 5.14 5.52 6
15.4% (4), 29.1% (2) 27.3% (6) + 6.4% (6) + 9.2% (6)
38% 45% + 17.6% + 23.9%
C2h D3d D3d D3d
DFT SW EDIP T3
10.12 10.8 10.2 11.6
27.7% (6) 30.9% (6) + 8.6% (6) + 11.4% (6)
27.3% 33.5% + 13.5% + 18.2%
D3d D3d D3d D3d
structure are listed in Table 1. Our DFT results 15.4% (4) and 29.1% (2) are reasonably in agreement with other LP calculations [ 13.8% (4), 29.03% (2) (Ref. [3]) and 12.2% (4), 25.8% (2) (Ref. [5])]. The slight differences may be induced by the use of shorter a0 of 5.48 A˚ [3] and ˚ In the three empirical potential 5.43 A˚ [5] than ours (5.494 A). calculations, the six nearest-neighbor atoms undergo a breathing model distortion, where the atoms move along the bonds between the nearest-neighbor atoms and their corresponding vacancy sites. Thus, the D3d symmetry of the initial states is not changed in the final states. In addition, SW leads to an inward distortion, which brings the nearest-neighbor atoms of each vacant site closer to each other to 2.826 A˚ (27.3%) and reduces the volume by 45%. Unlike SW, EDIP and T3 give an outward distortion, which makes the distances 6.4% and 9.2% longer than the unrelaxed structures. In addition, their volumes increase by 17.6% and 23.9%, respectively. The outward distortion has also been observed by Balamane et al. [20] in the V1 calculations. Although some prior Green’s function calculations [36,37], which lead to an outward relaxation of V1, Probert and Payne [35] found that the spurious outward relaxation was caused by the very small basis set size used in the calculations. The formation energies of relaxed V2 (Ef2) calculated by DFT, SW, EDIP and T3 are also listed in Table 1. Our DFT result 5.59 eV is approximately equal to the previous TB calculation 5.68 eV [38]. Wixom and Wright [4] have performed a series of DFT/GGA calculations in the supercells ranging from 216 to 1000 atoms and extra˚ Their result polated Ef2 for infinite-sized supercell at a0 ¼5.469 A. 5.363 eV is 0.227 eV less than our work. In addition, the binding energy of Wixom and Wright [4] (1.85 eV) is 0.06 eV larger than our result 1.79 eV. Our value is closer to the experimental estimate of Z1.6 eV [1]. From Table 1, it can be noted that, EDIP gives a perfect performance in the Ef2 calculation, which differs by 0.07 eV with our DFT result. In addition, the results of SW and EDIP are 0.45 less and 0.41 larger than our DFT value, respectively. In view of the accuracy of the empirical potentials, these two differences are reasonably admissible. The results shown in Table 1 indicate that empirical potentials give an invalid description of the symmetries in comparison with DFT. For the three potentials, only SW gives an inward distortion around V2 and its relaxation volume is in agreement with our DFT result. Moreover, its distortion magnitude is also close to the averaged DFT value of 20%. For EDIP and T3, the outward relaxations make the nearest-neighbor atoms away from each other. On the other hand, all the three potentials are competent for the calculation of Ef2. However, in consideration with structural properties, only the SW can be applied for the V2 calculations except for the case of symmetry.
3.2. Hexavacancy V6 is formed by removing six Si atoms from an hexagonal ring in the perfect crystal, which is illustrated in Fig. 2(b). In the initial state, each vacant site has two nearest-neighbors. It is consistent with earlier first-principle studies [12–14], the two nearest-neighbor atoms of each vacant site pair two by two in the D3d symmetry. In the
optimized geometry, the distances between these paired atoms are ˚ and the volume shortened to 2.808 A˚ (the ideal distance is 3.885 A) of V6 decreases by 27.3%. Our distortion magnitude is 0.052 A˚ ˚ of Akiyama and smaller than the DFT/GGA result 2.86 A˚ (3.88 A) Oshiyama [13]. Staab et al. [15] have performed the DFTB calcula˚ In the SW, EDIP and T3 tions and gave the similar value 2.7 A. calculations, all of them give a breathing model distortion. This radial deformation around the vacant sites leads to the D3d symmetry of the final states. For the SW, the relative distortion magnitude is 30.9% and the relative relaxation volume is 33.5%, which are similar to our DFT result of 27.7% and 27.3%, respectively. But for EDIP and T3, both of them undergo an outward relaxation, which are contrary to our DFT result. The formation energies of relaxed V6 (Ef6) calculated by DFT, SW, EDIP and T3 are listed in Table 2. In the earlier DFTB calculation, Staab et al. [15] gave the Ef6 of 10.5 eV, which are 0.38 eV larger than our DFT result 10.12 eV. Makhov and Lewis [12] have calculated Ef6 using LDA in the 216 supercell. Their result 9.41 eV is 0.71 eV lower than our GGA value. To compare the present work with the DBC model [11], which did not consider the relaxation energy, we calculated the formation energy of unrelaxed structure using the DFT method. The result of this work 13.9 eV agrees with the CDB value of 14.1 eV [11]. From Table 2, it can be seen that Ef6 of EDIP and SW is in agreement with DFT, where the differences limit to 6.7%. For T3, the formation energy 11.4 eV is larger than our DFT value by about 12.6%. This implies that the error of T3 result is greater than SW and EDIP. Comparison of the results of V6 in Table 2 indicates that SW does a good job both in formation energy and structural properties (distortion magnitude, relaxation volume and symmetry). For EDIP, although it provides a faithful description of Ef6, the outward relaxation restricts its application. In addition, T3 is the worst in the three classical potentials, which gives a bigger error in energy formation and fails in distortion magnitude. Therefore, it can be concluded that SW seems to be the best potential for the V6 calculations.
4. Summary In summary, the DFT method may provide accurate descriptions of vacancies in comparison with experiments and other firstprinciples calculations. For the empirical potentials, all of them have poor transferability in the JT distortion and undergo a breathing model distortion during the relaxations. For this reason, the potentials are failed in another quantum-mechanical phenomenon: the buckled dimer reconstruction on the Si(0 0 1) surface [20,39]. These unsatisfactory behaviors should be directly related to the form of potential. The short-range empirical potentials are fitted to experimental or DFT results, which are usually used to reproduce physical mechanisms or configurations at equilibrium. Since vacancies of Si involve rebonding and large atomic deformation, their atomic structures are far from the corresponding equilibrium states. This will probably require the atomic interactions more than three-body in the potential. Due to the poor transferability in the JT distortion, the three classical potentials give
470
C. Wang et al. / Physica B 406 (2011) 467–470
incorrect symmetry results for relaxed V2. In addition, EDIP and T3 give an outward relaxation in all cases, which is contrary to our DFT and SW directions. All the above results indicate that EDIP and T3 cannot be applied to describe the structural properties. On the other hand, although EDIP and T3 do good jobs in Ef2, SW should be the best potential V2 calculations since it gives good results for distortion magnitudes, relaxation volumes and formation energies. In the V6 calculations, all the three empirical potentials give the same symmetry of relaxed states as our DFT result. Even if both SW and EDIP provide a faithful description of Ef6, the outward relaxation of EDIP and T3 indicate that SW is the best potential. The main purpose of this work is the determination of the best empirical potentials to describe vacancies. In general, each potential has its own advantages and limitations for modeling vacancies. Despite their disadvantages, SW is relatively a better classical potential to describe V2 and V6. This may be explained by the cut-off radius of each potential. As stated in the above paragraph, the short-range potentials have limitations during investigation of ˚ of SW may vacancies. The relatively larger cut-off radius (3.77 A) ˚ and T3 result in good performance compared with EDIP (3.12 A) ˚ The EDIP and T3 have similar cut-off function expressions (3.0 A). with relatively small cut-off radii, which mean that they deal with fewer neighbors than SW in the simulation. This may cause that they are inapplicable in the structure calculations of V2 and V6.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
Acknowledgements This work was supported by the Fundamental Research Funds for the Central Universities and the NSF of China under Grant no. 10772062. References [1] G.D. Watkins, J.W. Corbett, Phys. Rev. 138 (1965) A543. [2] Y. Nagai, K. Inoue, Z. Tang, I. Yonenaga, T. Chiba, M. Saito, M. Hasegawa, Physica B 340–342 (2003) 518.
[29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
J.-I. Iwata, K. Shiraishi, A. Oshiyama, Phys. Rev. B 77 (2008) 115208. R.R. Wixom, A.F. Wright, Phys. Rev. B 74 (2006) 205208. S. Ogut, J.R. Chelikowsky, Phys. Rev. B 64 (2001) 245206. S. Ogut, J.R. Chelikowsky, Phys. Rev. Lett. 83 (1999) 3852. S.K. Estreiher, Phys. Status Solidi (b) 217 (2000) 513. J.P. Hirth, J. Lothe, Theory of Dislocations, Krieger, Malabar, 1992. B.L. Sopori, K. Jones, X.J. Deng, Appl. Phys. Lett. 61 (1992) 2560. B. Hourahine, R. Jones, A.N. Safonov, S. Oberg, P.R. Briddon, S.K. Estreicher, Phys. Rev. B 61 (2000) 12594. D.J. Chadi, K.J. Chang, Phys. Rev. B 38 (1988) 1523. D.V. Makhov, L.J. Lewis, Phys. Rev. Lett. 92 (2004) 255504. T. Akiyama, A. Oshiyama, J. Phys. Soc. Jpn. 70 (2001) 1627. J.L. Hastings, S.K. Estreicher, P.A. Fedders, Phys. Rev. B 56 (1997) 10215. T.E.M. Staab, A. Sieck, M. Haugk, M.J. Puska, Th. Frauenheim, H.S. Leipner, Phys. Rev. B 65 (2002) 115210. C.X. Li, Q.Y. Meng, K.Y. Zhong, C.Y. wang, Phys. Rev. B 77 (2008) 045211. C.Y. wang, Q.Y. Meng, Phys. Status Solidi RRL 3 (2009) 73. J.F. Justo, M.Z. Bazant, E. Kaxiras, V.V. Bulatov, S. Yip, Phys. Rev. B 58 (1998) 2539. M.S. Duesbery, B. Joos, D.J. Michel, Phys. Rev. B 43 (1991) 5143. H. Balamane, T. Halicioglu, W.A. Tiller, Phys. Rev. B 46 (1992) 2250. T.W. Poon, S. Yip, P.S. Ho, F.F. Abraham, Phys. Rev. Lett. 65 (1990) 2161. L.J. Porter, S. Yip, M. Yamaguchi, H. Kaburaki, M. Tang, J. Appl. Phys. 81 (1997) 96. J. Godet, L. Pizzagalli, S. Brochard, P. Beauchamp, J. Phys.: Condens. Matter 15 (2003) 6943. F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. M.Z. Bazant, E. Kaxiras, J.F. Justo, Phys. Rev. B 56 (1997) 8542. J. Tersoff, Phys. Rev. B 39 (1989) 5566. D. Sanchez-Portal, P. Ordejon, E. Artacho, J.M. Soler, Int. J. Quantum Chem. 65 (1997) 453. J.M. Soler, E. Artacho, J.D. Gale, A. Garcı´a, J. Junquera, P. Ordejo´n, D. Sa´nchezPortal, J. Phys.: Condens. Matter 14 (2002) 2745. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. L. Kleinman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, second ed., Cambridge University Press, Cambridge, 1992. C. Kittel, Introduction to Solid State Physics, Wiley, New York, 1986. M.I.J. Probert, M.C. Payne, Phys. Rev. B 67 (2003) 075204. O. Gunnarsson, O. Jepsen, O.K. Andersen, Phys. Rev. B 27 (1983) 7144. M. Scheffler, J.P. Vigneron, G.B. Bachelet, Phys. Rev. B 31 (1985) 6541. E.G. Song, E. Kim, Y.H. Lee, Y.G. Hwang, Phys. Rev. B 48 (1993) 1486. F. Liu, M.G. Lagally, Phys. Rev. Lett. 76 (1996) 3156.
Contents lists available at ScienceDirect
Physica B journal homepage: www.elsevier.com/locate/physb
Comparative study of the empirical interatomic potentials and density-functional simulations of divacancy and hexavacancy in silicon Chaoying Wang a,n, Zhenqing Wang a, Qingyuan Meng b a b
College of Aerospace and Civil Engineering, Harbin Engineering University, Harbin 150001, PR China Department of Astronautical Science & Mechanics, Harbin Institute of Technology, Harbin 150001, PR China
a r t i c l e in f o
abstract
Article history: Received 8 April 2010 Received in revised form 4 October 2010 Accepted 5 November 2010
The comparison between three classical potentials and density functional theory (DFT) is performed mainly for divacancy and hexavacancy in Si crystal. According to their performances on the formation energies and structural properties (distortion magnitude, relaxation volume and symmetry), the limitations and validities of classical potentials are discussed. It is found that the outward relaxation directions of EDIP and Tersoff (T3) are contrary to the DFT and Stillinger–Weber (SW) directions (inward), which restrict their application in the structural property calculations. Except for the divacancy symmetries, the results of SW are in agreement with the DFT results. Thus, it can be concluded that SW should be the best potential to describe V2 and V6. & 2010 Elsevier B.V. All rights reserved.
Keywords: Vacancies DFT Empirical potential
1. Introduction Vacancies, such as divacancy (V2) and ring hexavacancy (V6), are fundamentally and technologically important defects in silicon. V2 is particularly attractive from the experimental point of view because it can be easily produced by electron irradiation, and moreover, it is quite stable. Nowadays, the structural properties or formation energies of V2 with a variety of charge states have been explored in the experiments [1,2] and ab initio calculations [3–6]. V6 is the smallest ring vacancy cluster. It is likely to be an efficient gettering center for a range of impurities, such as carbon, oxygen and copper atoms [7]. In addition, it can also be the precursors of extended defects, such as prismatic dislocations [8], H-related platelets [9] or photoluminescent (PL) center [10]. In the pioneer work, Chadi and Chang [11] used the dangling bonds counting (DBC) model and predicted that V6 should be more stable than other types of vacancy clusters. This was verified in the subsequent density functional theory (DFT) [12–14] and density functional based tight binding (DFTB) [15] calculations. This type of simple model of V6 is important since it can give a rough idea about how energetically favorable structures may look like. Moreover, the most important vacancy aggregates are those that are the most stable, as they are the most likely ones to survive high-temperature treatments. On the other hand, it is noted that the theoretical interest of vacancies is mainly focused on the DFT calculations. However, this method is unable to deal with the large-scale vacancies-related
problems, such as the interactions of vacancies with extended defects (grain boundary, dislocation [16,17] and stacking fault), self-interstitial migration in large systems and so on, due to large atomic system and long time period calculations. A possible alternative for these cases might be the use of empirical potentials, which are computationally much less expensive. However, the accuracy of those potentials cannot be ensured. Thus, for a particular defect, it is necessary to test the validity before their applications. For this purpose, a lot of tests has been done previously, e.g., dislocations [18,19], elastic constants [18,20], surfaces [20,21], thermo-mechanical properties [22] and so on. Moreover, Godet et al. [23] compared the Stillinger–Weber (SW) [24], environment-dependent interatomic potential (EDIP) [25] and Tersoff (T3) [26] with ab initio method under large shear. For vacancies, monovacancy (V1) is widely used for testing the ability of empirical potentials [18,20]. But for V2 and V6, to our knowledge, little work has been done on the determination of the best interatomic potentials. In this work, the structural properties and formation energies of V2 and V6 are calculated by the use of DFT and other three classical potentials. The limitations and validities of empirical potentials for modeling vacancies are discussed by means of the comparison of the calculated results. For the comparison between the DFT and empirical potentials results, only the neutral vacancies are considered in the present work.
2. Computational methods n
Corresponding author. Tel.: + 86 045182519956. E-mail address: [email protected] (C. Wang).
0921-4526/$ - see front matter & 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.physb.2010.11.014
The first-principle calculations of the vacancies are conducted using SIESTA code [27,28]. The Perdew, Burke and Ernzerhof (PBE)
468
C. Wang et al. / Physica B 406 (2011) 467–470
[29] formulation of the generalized gradient approximation (GGA) is used for the electronic exchange and correlation. Core electrons are replaced by the non-local norm-conserving Troullier–Martins pseudopotential [30], which is factorized in the Kleinman– Bylander form [31]. Due to the requirement of more precise results, a real-space grid equivalent to 160 Ry plane wave cut-off is used. For the Brillouin zone (BZ) sampling, Monkhorst–Pack (MP) [32] mesh q3 shifted by (k0, k0, k0) is utilized. In all cases, atomic geometries are relaxed using the conjugate-gradient (CG) [33] ˚ method until interatomic forces are less than 0.01 eV/A. During the DFT calculations, the single-zeta (SZ), double-zeta (DZ) and doubled-zeta with polarization (DZP) basis sets can be used for Si atoms. The lattice constants (a0) with different basis sets ˚ respectively, which are are converged to 5.643, 5.556 and 5.494 A, 3.8%, 2.3% and 1.2%, within the experimental value 5.43 A˚ [34]. To remove spurious electronic interactions, the periodic supercell method is employed in the DFT calculations. With the different supercell sizes and a0 obtained, the total energies per atom with respect to q are plotted in Fig. 1. For all cases, q varies from 1 to 4, which corresponds to k0 ¼0 for q¼1 and k0 ¼0.5 for q¼2, 3, 4. Fig. 1 shows that the energy converges to different values with corresponding basis set. This implies that, except for computational techniques adopted, different a0 may lead to diverse results. In this work, the energy converges to 0.9 meV (DZ) and 0.04 meV (DZP) at q¼2 for 216 supercell. In addition, the mesh of 23 corresponds to a density of 0.031 A˚ 1, which is sufficient to make the total energy converge eventually [35]. In view of previous successful DFT calculations of V2 [4] and V6 [12,13], the 216 atoms supercell can be well converged eventually. Thus, the DFT calculations using the 216 supercell and DZP basis set at q¼2 are conducted in this work. Because of the great technological importance of Si, there are more than 30 empirical potentials, which have been developed and applied to different situations. In consideration to the advantages and limitations of each potential, the widely used and tested SW and T3 potentials [24] were applied for this work. Moreover, EDIP is known to describe both 301- and 901-partial dislocations in the glide set {1 1 1} [18], which was also used in this work. In fact, the SW potential is fitted to experimental properties of the diamondcubic (dc) and molten phases of silicon. Its formulation consists of a linear combination of two- and three-body terms and contains eight parameters. The EDIP also includes two- and three-body terms. which depend on the local atomic environment given by an effective coordination number. Its thirteen parameters in the functional form are obtained by fitting to a set of ab initio results,
Fig. 1. Convergence tests for the energy per atom as a function of Brillouin zone sampling for different sized supercells and basis sets.
which include various bulk phases and defect structures. The Tersoff potential consists of many-body interactions included in a bond order term. In addition, the eleven parameters are fitted to ab initio results for several Si polytypes. It has three versions in the history and the last version (T3) is used in this work. The total energies are calculated by the molecular dynamics (MD) method with the NPT (the atoms number, pressure and temperature are kept constant) ensemble and a time step of 1 fs. The vacancy is formed in the middle of MD model, which is 9.9 nm along the three vectors [1 0 0], [0 1 0] and [0 0 1]. In the calculations, the periodic boundary conditions (PBC) are applied in the three directions of the model. The total energy is minimized during the annealing process. First, atomic structure is relaxed for 10,000 steps at 0 K. And then, the temperature becomes 10 K higher at every 5000 steps until a desired temperature is achieved. By the reverse process, the temperature drops to 0 K and the energy is obtained. In the present work, zero-temperature formation energies of vacancies are computed using the expression Efn ¼ EðNnÞ
Nn EðNÞ N
ð1Þ
where E(N) is the total energy of the perfect crystal containing N atoms, E(N n) is the total energy of the crystal containing n vacant sites. As a measure of the variation in vacancy volumes resulting from the relaxation, we also calculate the relative relaxation volume
DO ¼ ðOO0 Þ=O0 100%
ð2Þ
where O0 and O are the volumes of the vacancy before and after relaxation, respectively.
3. Results and discussion 3.1. Divacancy Removing two adjacent atoms in the perfect crystal generates an unrelaxed V2, where atomic structures are of D3d symmetry. From Fig. 2(a), we can see that there are two sets of three distances between the three nearest-neighbors of each vacant site. During the DFT relaxation process, the Jahn–Teller (JT) distortion reduces the ideal D3d symmetry to C2h. Moreover, the pairing of the atoms 2–3 and 4–5 forms the large pairing (LP) [1] structure, where d2,3 ¼ d4,5 od1,2 ¼d1,3 ¼d4,6 ¼d5,6. During this process, the volume of V2 reduces by 38%. Except for the LP structure, Iwata et al. [3] claimed that the resonant bond (RB) structure was more stable than the DFT calculations in the supercell ranging from 64 to 1000 atoms. The RB structure is signaled by d2,3 ¼ d4,5 4d1,2 ¼d1,3 ¼d4,6 ¼d5,6. But in another large scale (1000-atom) DFT calculation [4], Wixom and Wright found that the LP structure was more stable than the RB. This conclusion was proved by other DFT calculations [5,6] and experiment [2]. The six distortion amplitudes relative to the unrelaxed
Fig. 2. Schematic view of (a) divacancy (V2), (b) hexavacancy (V6) and their nearestneighbor atoms.
C. Wang et al. / Physica B 406 (2011) 467–470
Table 1 Structural properties and formation energies (Ef2) of relaxed divacancy. Dd (%) indicates the fractional distortion magnitudes, which is relative to the ideal distance ˚ The signs and + indicate the inward and outward relaxation, respectively. 3.885 A.
469
Table 2 Structural properties and formation energies (Ef6) of relaxed hexavacancy. Dd (%) indicates the fractional distortion magnitudes, which is relative to the ideal distance ˚ The signs and + indicate the inward and outward relaxation, respectively. 3.885 A.
Method
Ef2 (eV)
Dd (%)
DO (%)
Symmetry
Method
Ef6 (eV)
Dd (%)
DO (%)
Symmetry
DFT SW EDIP T3
5.59 5.14 5.52 6
15.4% (4), 29.1% (2) 27.3% (6) + 6.4% (6) + 9.2% (6)
38% 45% + 17.6% + 23.9%
C2h D3d D3d D3d
DFT SW EDIP T3
10.12 10.8 10.2 11.6
27.7% (6) 30.9% (6) + 8.6% (6) + 11.4% (6)
27.3% 33.5% + 13.5% + 18.2%
D3d D3d D3d D3d
structure are listed in Table 1. Our DFT results 15.4% (4) and 29.1% (2) are reasonably in agreement with other LP calculations [ 13.8% (4), 29.03% (2) (Ref. [3]) and 12.2% (4), 25.8% (2) (Ref. [5])]. The slight differences may be induced by the use of shorter a0 of 5.48 A˚ [3] and ˚ In the three empirical potential 5.43 A˚ [5] than ours (5.494 A). calculations, the six nearest-neighbor atoms undergo a breathing model distortion, where the atoms move along the bonds between the nearest-neighbor atoms and their corresponding vacancy sites. Thus, the D3d symmetry of the initial states is not changed in the final states. In addition, SW leads to an inward distortion, which brings the nearest-neighbor atoms of each vacant site closer to each other to 2.826 A˚ (27.3%) and reduces the volume by 45%. Unlike SW, EDIP and T3 give an outward distortion, which makes the distances 6.4% and 9.2% longer than the unrelaxed structures. In addition, their volumes increase by 17.6% and 23.9%, respectively. The outward distortion has also been observed by Balamane et al. [20] in the V1 calculations. Although some prior Green’s function calculations [36,37], which lead to an outward relaxation of V1, Probert and Payne [35] found that the spurious outward relaxation was caused by the very small basis set size used in the calculations. The formation energies of relaxed V2 (Ef2) calculated by DFT, SW, EDIP and T3 are also listed in Table 1. Our DFT result 5.59 eV is approximately equal to the previous TB calculation 5.68 eV [38]. Wixom and Wright [4] have performed a series of DFT/GGA calculations in the supercells ranging from 216 to 1000 atoms and extra˚ Their result polated Ef2 for infinite-sized supercell at a0 ¼5.469 A. 5.363 eV is 0.227 eV less than our work. In addition, the binding energy of Wixom and Wright [4] (1.85 eV) is 0.06 eV larger than our result 1.79 eV. Our value is closer to the experimental estimate of Z1.6 eV [1]. From Table 1, it can be noted that, EDIP gives a perfect performance in the Ef2 calculation, which differs by 0.07 eV with our DFT result. In addition, the results of SW and EDIP are 0.45 less and 0.41 larger than our DFT value, respectively. In view of the accuracy of the empirical potentials, these two differences are reasonably admissible. The results shown in Table 1 indicate that empirical potentials give an invalid description of the symmetries in comparison with DFT. For the three potentials, only SW gives an inward distortion around V2 and its relaxation volume is in agreement with our DFT result. Moreover, its distortion magnitude is also close to the averaged DFT value of 20%. For EDIP and T3, the outward relaxations make the nearest-neighbor atoms away from each other. On the other hand, all the three potentials are competent for the calculation of Ef2. However, in consideration with structural properties, only the SW can be applied for the V2 calculations except for the case of symmetry.
3.2. Hexavacancy V6 is formed by removing six Si atoms from an hexagonal ring in the perfect crystal, which is illustrated in Fig. 2(b). In the initial state, each vacant site has two nearest-neighbors. It is consistent with earlier first-principle studies [12–14], the two nearest-neighbor atoms of each vacant site pair two by two in the D3d symmetry. In the
optimized geometry, the distances between these paired atoms are ˚ and the volume shortened to 2.808 A˚ (the ideal distance is 3.885 A) of V6 decreases by 27.3%. Our distortion magnitude is 0.052 A˚ ˚ of Akiyama and smaller than the DFT/GGA result 2.86 A˚ (3.88 A) Oshiyama [13]. Staab et al. [15] have performed the DFTB calcula˚ In the SW, EDIP and T3 tions and gave the similar value 2.7 A. calculations, all of them give a breathing model distortion. This radial deformation around the vacant sites leads to the D3d symmetry of the final states. For the SW, the relative distortion magnitude is 30.9% and the relative relaxation volume is 33.5%, which are similar to our DFT result of 27.7% and 27.3%, respectively. But for EDIP and T3, both of them undergo an outward relaxation, which are contrary to our DFT result. The formation energies of relaxed V6 (Ef6) calculated by DFT, SW, EDIP and T3 are listed in Table 2. In the earlier DFTB calculation, Staab et al. [15] gave the Ef6 of 10.5 eV, which are 0.38 eV larger than our DFT result 10.12 eV. Makhov and Lewis [12] have calculated Ef6 using LDA in the 216 supercell. Their result 9.41 eV is 0.71 eV lower than our GGA value. To compare the present work with the DBC model [11], which did not consider the relaxation energy, we calculated the formation energy of unrelaxed structure using the DFT method. The result of this work 13.9 eV agrees with the CDB value of 14.1 eV [11]. From Table 2, it can be seen that Ef6 of EDIP and SW is in agreement with DFT, where the differences limit to 6.7%. For T3, the formation energy 11.4 eV is larger than our DFT value by about 12.6%. This implies that the error of T3 result is greater than SW and EDIP. Comparison of the results of V6 in Table 2 indicates that SW does a good job both in formation energy and structural properties (distortion magnitude, relaxation volume and symmetry). For EDIP, although it provides a faithful description of Ef6, the outward relaxation restricts its application. In addition, T3 is the worst in the three classical potentials, which gives a bigger error in energy formation and fails in distortion magnitude. Therefore, it can be concluded that SW seems to be the best potential for the V6 calculations.
4. Summary In summary, the DFT method may provide accurate descriptions of vacancies in comparison with experiments and other firstprinciples calculations. For the empirical potentials, all of them have poor transferability in the JT distortion and undergo a breathing model distortion during the relaxations. For this reason, the potentials are failed in another quantum-mechanical phenomenon: the buckled dimer reconstruction on the Si(0 0 1) surface [20,39]. These unsatisfactory behaviors should be directly related to the form of potential. The short-range empirical potentials are fitted to experimental or DFT results, which are usually used to reproduce physical mechanisms or configurations at equilibrium. Since vacancies of Si involve rebonding and large atomic deformation, their atomic structures are far from the corresponding equilibrium states. This will probably require the atomic interactions more than three-body in the potential. Due to the poor transferability in the JT distortion, the three classical potentials give
470
C. Wang et al. / Physica B 406 (2011) 467–470
incorrect symmetry results for relaxed V2. In addition, EDIP and T3 give an outward relaxation in all cases, which is contrary to our DFT and SW directions. All the above results indicate that EDIP and T3 cannot be applied to describe the structural properties. On the other hand, although EDIP and T3 do good jobs in Ef2, SW should be the best potential V2 calculations since it gives good results for distortion magnitudes, relaxation volumes and formation energies. In the V6 calculations, all the three empirical potentials give the same symmetry of relaxed states as our DFT result. Even if both SW and EDIP provide a faithful description of Ef6, the outward relaxation of EDIP and T3 indicate that SW is the best potential. The main purpose of this work is the determination of the best empirical potentials to describe vacancies. In general, each potential has its own advantages and limitations for modeling vacancies. Despite their disadvantages, SW is relatively a better classical potential to describe V2 and V6. This may be explained by the cut-off radius of each potential. As stated in the above paragraph, the short-range potentials have limitations during investigation of ˚ of SW may vacancies. The relatively larger cut-off radius (3.77 A) ˚ and T3 result in good performance compared with EDIP (3.12 A) ˚ The EDIP and T3 have similar cut-off function expressions (3.0 A). with relatively small cut-off radii, which mean that they deal with fewer neighbors than SW in the simulation. This may cause that they are inapplicable in the structure calculations of V2 and V6.
[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
Acknowledgements This work was supported by the Fundamental Research Funds for the Central Universities and the NSF of China under Grant no. 10772062. References [1] G.D. Watkins, J.W. Corbett, Phys. Rev. 138 (1965) A543. [2] Y. Nagai, K. Inoue, Z. Tang, I. Yonenaga, T. Chiba, M. Saito, M. Hasegawa, Physica B 340–342 (2003) 518.
[29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
J.-I. Iwata, K. Shiraishi, A. Oshiyama, Phys. Rev. B 77 (2008) 115208. R.R. Wixom, A.F. Wright, Phys. Rev. B 74 (2006) 205208. S. Ogut, J.R. Chelikowsky, Phys. Rev. B 64 (2001) 245206. S. Ogut, J.R. Chelikowsky, Phys. Rev. Lett. 83 (1999) 3852. S.K. Estreiher, Phys. Status Solidi (b) 217 (2000) 513. J.P. Hirth, J. Lothe, Theory of Dislocations, Krieger, Malabar, 1992. B.L. Sopori, K. Jones, X.J. Deng, Appl. Phys. Lett. 61 (1992) 2560. B. Hourahine, R. Jones, A.N. Safonov, S. Oberg, P.R. Briddon, S.K. Estreicher, Phys. Rev. B 61 (2000) 12594. D.J. Chadi, K.J. Chang, Phys. Rev. B 38 (1988) 1523. D.V. Makhov, L.J. Lewis, Phys. Rev. Lett. 92 (2004) 255504. T. Akiyama, A. Oshiyama, J. Phys. Soc. Jpn. 70 (2001) 1627. J.L. Hastings, S.K. Estreicher, P.A. Fedders, Phys. Rev. B 56 (1997) 10215. T.E.M. Staab, A. Sieck, M. Haugk, M.J. Puska, Th. Frauenheim, H.S. Leipner, Phys. Rev. B 65 (2002) 115210. C.X. Li, Q.Y. Meng, K.Y. Zhong, C.Y. wang, Phys. Rev. B 77 (2008) 045211. C.Y. wang, Q.Y. Meng, Phys. Status Solidi RRL 3 (2009) 73. J.F. Justo, M.Z. Bazant, E. Kaxiras, V.V. Bulatov, S. Yip, Phys. Rev. B 58 (1998) 2539. M.S. Duesbery, B. Joos, D.J. Michel, Phys. Rev. B 43 (1991) 5143. H. Balamane, T. Halicioglu, W.A. Tiller, Phys. Rev. B 46 (1992) 2250. T.W. Poon, S. Yip, P.S. Ho, F.F. Abraham, Phys. Rev. Lett. 65 (1990) 2161. L.J. Porter, S. Yip, M. Yamaguchi, H. Kaburaki, M. Tang, J. Appl. Phys. 81 (1997) 96. J. Godet, L. Pizzagalli, S. Brochard, P. Beauchamp, J. Phys.: Condens. Matter 15 (2003) 6943. F.H. Stillinger, T.A. Weber, Phys. Rev. B 31 (1985) 5262. M.Z. Bazant, E. Kaxiras, J.F. Justo, Phys. Rev. B 56 (1997) 8542. J. Tersoff, Phys. Rev. B 39 (1989) 5566. D. Sanchez-Portal, P. Ordejon, E. Artacho, J.M. Soler, Int. J. Quantum Chem. 65 (1997) 453. J.M. Soler, E. Artacho, J.D. Gale, A. Garcı´a, J. Junquera, P. Ordejo´n, D. Sa´nchezPortal, J. Phys.: Condens. Matter 14 (2002) 2745. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865. N. Troullier, J.L. Martins, Phys. Rev. B 43 (1991) 1993. L. Kleinman, D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. H.J. Monkhorst, J.D. Pack, Phys. Rev. B 13 (1976) 5188. W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, second ed., Cambridge University Press, Cambridge, 1992. C. Kittel, Introduction to Solid State Physics, Wiley, New York, 1986. M.I.J. Probert, M.C. Payne, Phys. Rev. B 67 (2003) 075204. O. Gunnarsson, O. Jepsen, O.K. Andersen, Phys. Rev. B 27 (1983) 7144. M. Scheffler, J.P. Vigneron, G.B. Bachelet, Phys. Rev. B 31 (1985) 6541. E.G. Song, E. Kim, Y.H. Lee, Y.G. Hwang, Phys. Rev. B 48 (1993) 1486. F. Liu, M.G. Lagally, Phys. Rev. Lett. 76 (1996) 3156.