Tuning Coupling Behavior of Stacked Heterostructures Based on MoS2, WS2, and WSe2

The interlayer interaction of vertically stacked heterojunctions is very sensitive to the interlayer spacing, which will affect the coupling between the monolayers and allow band structure modulation. Here, with the aid of density functional theory (DFT) calculations, an interesting phenomenon is found that MoS2-WS2, MoS2-WSe2, and WS2-WSe2 heterostructures turn into direct-gap semiconductors from indirect-gap semiconductors with increasing the interlayer space. Moreover, the electronic structure changing process with interlayer spacing of MoS2-WS2, MoS2-WSe2, and WS2-WSe2 is different from each other. With the help of variable-temperature spectral experiment, different electronic transition properties of MoS2-WS2, MoS2-WSe2, and WS2-WSe2 have been demonstrated. The transition transformation from indirect to direct can be only observed in the MoS2-WS2 heterostructure, as the valence band maximum (VBM) at the Γ point in the MoS2-WSe2 and WS2-WSe2 heterostructure is less sensitive to the interlayer spacing than those from the MoS2-WS2 heterostructure. The present work highlights the significance of the temperature tuning in interlayer coupling and advance the research of MoS2-WS2, MoS2-WSe2, and WS2-WSe2 based device applications.

Semiconductor heterojunctions have played a significant role in the rich collection of unimagined electronic structures and optical properties 1 . The emergence of atom-thin transition-metal dichalcogenides (TMDs) as a new class of two dimensional semiconducting materials, which are almost as thin, transparent, and flexible as graphene with a wealth of new physical phenomena [2][3][4][5] , creates exciting new opportunities to push semiconductor heterostructures toward a new frontier 6,7 . Vertically stacked van der Waals TMDs heterostructures have been recognized as a powerful platform to create atomically thin heterostructures [8][9][10] . These heterojunctions have an optically active band gap with bound electrons and holes localized in individual monolayers 11,12 , and their energy and luminescence intensity are highly tunable with the applied vertical gate voltage, laser intensity, and annealing time 12,13 . A more desirable way to control the interlayer coupling properties is adjusting the interlayer spacing, which is closely related with the temperature 14,15 . Temperature induces lattice parameter changes and modulates the band structure, altering the optical properties of the heterojunctions. Further more, it is important to investigate the temperature induced phenomena in order to distinguish the phenomena influenced by other factors, such as the surface quality, the strain applied on the heterostructures, and even the number of atomic layers, etc. However, such an interesting and important research of variable-temperature tuning for TMDs heterostructures has not been reported up to now.
In this article, high-quality MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures have been fabricated by the polystyrene film transferred technique. They are systematically studied by tuning interlayer coupling with temperature for the first time. Furthermore, it is discovered that the temperature is an extremely sensitive factor to the interlayer coupling in comparison with the annealing tuning. Based on the density functional theory (DFT) calculations, an interesting phenomenon is found that MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures turn into direct-gap semiconductors from indirect-gap semiconductors with increasing the interlayer space. However, in the MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures, the change of the bandgap properties with interlayer spacing differs from each other. The results has been observed in the variable-temperature experiment. Our results not only compare the MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures tuning character, but also open up a new direction for 2D applications where external modulation of bandgap and optical properties is desired.

Results and discussion
Raman scattering. The vertically stacked heterostructures have been investigated by Raman scattering and PL spectroscopy. The optical images and Raman spectra are shown in Fig. 1 Fig. 1(e)] seem to be a simple superposition of the monolayer MoS 2 and WS 2 . Even after the annealing at 100 °C in argon gas for 6 hours, except for increasing the intensity [red line in Fig. 1 13 . The frequency difference between the E g 2 1 and A 1g Raman modes is 18.6 cm −1 for monolayer MoS 2 in MoS 2 -WSe 2 [green line in Fig. 1(j)], and 4.7 cm −1 for monolayer WSe 2 in WS 2 -WSe 2 [green line in Fig. 1(o)], respectively 19 . While in the WSe 2 -based heterostructures, an additional weak Raman peak pointed to interlayer coupling near 300 cm −1 can be observed [insert figures in Fig. 1(j) and (o)], which disappears in the monolayer WSe 2 . It corresponds to the B 2g resonance mode of WSe 2 20 . In general, the B 2g signature mode is only active in the bilayer or few-layer WSe 2 , which could reflect the presence of the additional interlayer interaction in the present TMDs heterostructures 20 . It demonstrates that the B 2g mode is more sensitive to interlayer interaction than the out-of-plane modes A 1g [19][20][21] .
Band structures for heterostructures. The conduction band minimum (CBM) and valence band maximum (VBM) at the K point are primarily composed of the Mo or W d z 2 states, while the CBM at K-Γ and the VBM at Γ are dominated by the Mo or W d z 2 and and S p states. The Mo d xz yz , states are far away from the Fermi energy 13 . Because of different orbital character, the highest valence and lowest conduction states respond very differently to the interlayer spacing and coupling. For a detailed research for the interlayer spacing and coupling, the DFT calculations for MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures have been carried out, which are shown in Fig. 2. The equilibrium interlayer separation distances with geometry optimization for MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures are 6.155 Å, 6.418 Å, and 6.702 Å, respectively. The corresponding band structures near the equilibrium layer spacing are shown in Fig. 2 . Interestingly, at the equilibrium layer spacing of 6.155 Å in the MoS 2 -WS 2 heterostructure, the VBM is at the Γ point because of the interaction between layers [ Fig. 2(b)]. This results in an indirect transition from the VBM at the Γ point to the CBM at the K points. It is different from the direct transition situated at K point of the monolayer TMDs. The VBM at the Γ point becomes even higher to the K point when the interlayer separation distances is reduced to 5.655 Å in Fig. 2(a). The results are more obvious in Fig. 2(d), which shows the variations of the highest valence and lowest conduction states as functions of the layer spacing in MoS 2 -WS 2 . It can be seen that the reduction in the interlayer spacing modifies the VBM and CVM values at high symmetry points effectively. As the VBM at the Γ point involving p-d orbital coupling, it changes dramatically. As a result, the indirect bandgap increases obviously for shorter layer spacing. With decreasing the layer spacing from 7.155 Å to 5.155 Å, not only the VBM at Γ point becomes higher with respect to the K point, but also the CBM at a midpoint between K and Γ points becomes lower than that at the K point. However, the direct excitonic transition energy at the K point shows small changes.
The  the layer spacing is (a) 0.5 Å smaller than the equilibrium distances, (b) equal to the equilibrium distances, (c) and 0.5 Å more than the equilibrium distances. (e-g) Band structures for MoS 2 -WSe 2 when the layer spacing is (e) 0.5 Å smaller than the equilibrium distances, (f) equal to the equilibrium distances, and (g) 0.5 Å more than the equilibrium distances. (i-k) Band structures for WS 2 -WSe 2 when the layer spacing is (i) 0.5 Å smaller than the equilibrium distances, (j) equal to the equilibrium distances, and (k) 0.5 Å more than the equilibrium distances. (d,h,l) Variation of the VBM and CBM with respect to the the layer spacing for (d) MoS 2 -WS 2 , (h) MoS 2 -WSe 2 , and (l) WS 2 -WSe 2 . VBM K and CBM K refer to the local highest valence and lowest conduction states at K point, VBM Γ to the Γ point, and CBM K−Γ to the path K-Γ . only 1 eV, which is less than those of about 1.3 eV in the MoS 2 -WS 2 heterostructure. This phenomenon indicates that the interlayer coupling tuning with interlayer separation distances for WS 2 -WSe 2 heterostructure is relatively difficult, which is in accordance with the presenting experimental results.
Variable-temperature tuning process. From tight-binding theory or quantum tunneling model, the interlayer interaction is expected to be exponentially sensitive to the interlayer distance. In order to further confirm the theoretical calculation and investigate the interlayer coupling with different interlayer distance, variable-temperature tuning process were performed in MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures. It has been known that the physical properties of the heterostructures could be affected by the twist angle between the two layers as observed in graphene and graphitic materials 23,24 . To avoid the interference of the twist angle, the whole variable-temperature experiment was carried on the same point of the heterostructure. For the purpose of obtaining information on the band properties of the heterojunctions, optical properties are characterized with PL spectroscopies in Fig. 3.
All MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures form a type II band alignment 22,25 . In type II heterojunctions, the CBM and VBM reside in two separate materials. Moreover, the type II band alignment and built-in potential in the heterojunctions can facilitate the photoexcited electron-hole separation and lead to an enhanced photoswitching performance compared to that in MoS 2 WS 2 and WSe 2 26 . For example, in MoS 2 -WS 2 heterostructure, due to the type-II band alignment, photoexcited electrons and holes will relax (dashed lines in Fig. 3(b)) to the conduction band edge of MoS 2 and the valence band edge of WS 2 , respectively. The coulomb attraction between electrons in MoS 2 and holes in WS 2 gives rise to an interlayer exciton P hetero . It is analogous to spatially indirect excitons in coupled quantum wells 27 . Therefore, a new peak at 1.77 eV is observed P [ hetero in Fig. 3(a)], lying interestingly at a lower energy than the peak for the two constituent single layers. It is noteworthy that the relative intensity for P hetero and MoS 2 , WS 2 changed after the annealing progress. The relative intensity of P hetero to P MoS 2 is expressed as I h Mo / , P hetero to P W S 2 is I h W / . After the annealing progress, I h Mo / increases from 0.65 to 0.87, I h W / increases from 0.3 to 0.61. Such a pronounced luminescence effect in the heterostructure suggests that most electrons excited in WS 2 transfer to the lower states in MoS 2 , instead of forming excitons in WS 2 and recombine radiatively. The annealing progress promote the electrons transfer in the heterostructure, as the annealing is able to drive out trapped residual molecules 28 . However, in equilibrium condition of the MoS 2 -WSe 2 and WS 2 -WSe 2 heterostructures, the direct gaps at the K point are around 1 eV [ Fig. 2(f) and (j)].which are different from the indirect transition of MoS 2 -WS 2 22 . Under these conditions, the PL spectra from 1.4 eV to 2.3 eV for MoS 2 -WSe 2 and WS 2 -WSe 2 heterostructures in Fig. 3(c) and (d) appear to comprise the addition of the constituent layers.
In the variable-temperature tuning process for the MoS 2 -WS 2 , the excitonic transitions of MoS 2 (P MoS 2 ), WS 2 (P W S 2 ), and interlayer exciton (P hetero ) are fitted with the integrated intensity peak P M , with the purpose of describing the intensity change with temperature. Figure 4(a) shows the luminescence evolutions of the MoS 2 -WS 2 heterostructure as a function of temperature from 77 to 275 K with a stability of about 0.1 K. When the excitonic transitions of MoS 2 (P MoS 2 ), WS 2 (P W S 2 ), and interlayer exciton (P hetero ) are fitted with the integrated intensity peak P M , the PL spectra can be deconvoluted into three major peaks when the temperature is below 225 K. The resonance at about 2 eV corresponds to the B exciton, which attributes to the top of valence-band splitting due to the strong spin-orbital interaction in MoS 2 29 . In addition to the above excitons, another transition at a lower energy can be observed at about 1.7 eV, which seems different from the previous work 30,31 . The transition becomes negligible in the temperature above 250 K [P indiret in Fig. 4(a),(d) and (g)]. This characteristic is profoundly different from the behaviors of other double layer transition metal sulfides. Nanostructures obtained from indirect bandgap transition metal sulfides emit an indirect transition excitons upon photoexcitation, but the luminescence is present in the room temperature as well ( Figure S1) 32 . In order to identify the peak P indiret , another significant change should be noticed. When the temperature increases above 250 K, the intensity of the peak P M increases dramaticly by about three times. The integral PL intensity generally decreases with increasing the temperature. The phenomena can be observed in temperature dependent PL measurement of monolayer MoS 2 , as shown in Figure S2. This quenching effect is due to that the nonradiative electron-hole recombination rate increases exponentially with increasing the temperature 33,34 . The nonradiative channels, such as trapping by surface/defect/ ionized impurity states, become thermally activated with increasing the temperature and the nonradiative lifetime τ NR is expressed as ref. 35 where τ 0 is the pre-exponential factor, E a is the activation energy in the thermal quenching process and k B is the Boltzmann constant. When the temperature increases, the decreasing τ NR results in a decrease in the luminescence intensity, which can be expressed as ref. 36 where I 0 is the intensity at 0 K, A is a parameter related to radiative lifetime (τ R ) as A = τ R /τ 0 , and τ R is the radiative lifetime. However, the peak P M reveals the monotonic decrease of PL intensity with increasing the temperature only when the temperature is below 250 K, and the PL intensity is even stronger when the temperature is above 250 K in spite of the increasing temperature. In order to explain the extraordinary change of the P M intensity, the luminescence quantum efficiency is introduced. As the luminescence physics mechanism for transition metal sulfide is different from traditional silicon nanocrystals. In silicon nanocrystals, the photoluminescence originates from quantum confined electronic states with increased emission energy at decreased nanoparticle size. The radiative transition rate is quite low. In contrast, luminescence for transition metal sulfide arises from electronic transitions, which shows a much higher radiative recombination rate. Therefore, the extraordinary change of the P M intensity is explained by luminescence quantum efficiency, which is approximated by ref. 37 where k rad is the rates of radiative recombination, k relax is electron relaxation within the conduction and valence bands, and k defect is defect trapping, respectively. The k red is not likely to change appreciably with the change of the interlayer spacing in heterostructure, because the direct excitonic transitions at Γ point do not show significant change. The k defect is stable relatively. Therefore, when the temperature increases above 250 K, the significant enhanced luminescence of P M has to be attributed to the dramatically reduced electronic relaxation k relax . The decrease of interband relaxation rate strongly suggests a substantial change in the MoS 2 -WS 2 heterostructure electronic structure when the temperature increases from 225 K to 250 K. It is noteworthy that the value of the phonon-assisted k relax is very large for the indirect bandgap semiconductor. Because when the indirect bandgap disappears, the decay rate via phonons decreases enormously 32,37 . As a consequence, when the temperature increases from 225 K to 250 K, the indirect bandgap disappears. The substantial change is in accordance with the transition transformation analysed in the DFT calculations. The indirect transition from the VBM at Γ point to the CBM at K point disappeared with increasing the interlayer spacing. Based on the above analysis, the transition P indiret disappeared above 250 K is attributed to indirect band gap emission. It is noteworthy that, the MoS 2 -WSe 2 prefers direct transition at the equilibrium layer spacing, and the VBM at the Γ point for the MoS 2 -WSe 2 band structure changes less obviously than the MoS 2 -WS 2 [ Fig. 2(e-h)]. In the luminescence evolutions of MoS 2 -WSe 2 heterostructure [ Fig. 4(b)], another indirect transition peak P indiret emerges at the temperature below 225 K [ Fig. 4(e)]. However, the luminescence trends for the P MoS 2 and P WSe 2 are different from MoS 2 -WS 2 . Among the heterostructures studied, the atomic differences in MoS 2 -WSe 2 are the most obvious. Both the transition metals and the sulfur family elements are different. As a consequence, the WSe 2 attributed part in the band structure is significantly influenced in the heterostructure [blue line in Fig. 2(f)]. An indirect transition from the VBM of WSe 2 to the WSe 2 induced conduction band between Γ and K points is caused. Moreover, the P indiret is quenching with increasing the temperature [ Fig. 4b and e)], as the influence in the band structure of the WSe 2 attributed part declines when the interlayer spacing increases. Along with the quenching of the indirect transition, the direct transition excitons increase, and the P WSe 2 increases gradually. However, the VBM at the Γ point in the MoS 2 -WSe 2 heterostructure is less sensitive to the interlayer spacing than the MoS 2 -WS 2 heterostructure[ Fig. 2]. Therefore, the dramatically increasing intensity of the PL has not been observed here.
Not only for MoS 2 -WSe 2 , but the VBM at the Γ point for WS 2 -WSe 2 heterostructure is insensitive to the interlayer spacing [ Fig. 2(l)]. Moreover, the indirect transition for WSe 2 in the WS 2 -WSe 2 is less obvious than that in the MoS 2 -WSe 2 heterostructure. As a consequence, the P W S 2 and P WSe 2 present an redshift for conventional semiconductor. In addition, it is important to note that the B exciton peak of WS 2 and WSe 2 can be observed in the enhanced PL spectra of WS 2 -WSe 2 heterostructure in Fig. 4(i). The peaks near 1.98 eV and above 2.3 eV at 77 K are attributed to be the A and B excition of WS 2 17,32 . The peak near 2.13 eV at 77 K corresponds to the B excition of WSe 2 . It is 430 meV higher than the A excition, which is in agreement with the previous reports 32 . The B exciton peak for WS 2 and WSe 2 with such an intensity is rarely observed in the previous PL measurement. What is more, it is discovered for the first time that the B exciton peak for WSe 2 quenches with increasing the temperature.
Along with changing of the luminescent intensity, the direct transition in the heterostructures experiences a large redshift with increasing the temperature, as shown in Fig. 4(f-h). Such behavior is similar to the response of conventional semiconductors under high temperature, which result from the increased electron-phonon interactions and slight changes in bonding lengths. Thus, it provides a method to evaluate temperature of the semiconductor 38,39 . By employing a modified Varshni relationship, the temperature dependence of the PL peak position is fitted using 40,41 where E 0 is the emission energy at zero absolute temperature, S is the Huang-Rhys factor that represents the coupling strength of exciton-phonon, 〈 ħω〉 is the average phonon energy, ħ and k B are the Plancks and Boltzmann constant, respectively. As shown in Fig. 4(f-h), the fitting parameters are extracted. Table 1 lists the fitting values of E 0 , S, and 〈 ħω〉 for MoS 2 , WS 2 , and WSe 2 in the heterostructures. By comparing these parameters, especially S, between the TMD in the heterostructures and monolayer 38,42 , the contribution of the interlayer coupling can be informed. Moreover, with the help of this expression, it is possible to derive the temperature difference of the above heterostructure in the device where the PL is measured by comparing the emission energy difference 38,43 . However, from Fig. 4(g,h), the peak P indiret do not show an obvious redshift as the P M , P MoS 2 , or P WSe 2 , which is due to the gradual decreasing VBM at the Γ point caused by increasing the interlamellar spacing with the temperature. As shown in Figure S1, this decrease of VBN at the Γ point results in a blueshift of the indirect transition in double layer WS 2 .

Conclusions
To summarize, the interlayer coupling tuning of 2D heterostructures formed with CVD-grown monolayer MoS 2 , WS 2 , and WSe 2 is carried out by thermal annealing process and variable-temperature experiment. By comparing the tuning methods, the conclusion is made that the temperature is an extremely sensitive factor to the interlayer coupling. Based on the DFT calculations, an interesting phenomenon is found that MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures turn into direct-gap semiconductors from indirect-gap semiconductors with increasing the interlayer space. Furthermore in the MoS 2 -WS 2 , MoS 2 -WSe 2 , and WS 2 -WSe 2 heterostructures, the electronic structure changing process with interlayer spacing is different from each other. Our results highlight the significance of interlayer coupling in tuning the light emission of TMDs and offer a general route to prepare large-area TMD tandem structures for fundamental study as well as electronic and photovoltaic applications.

Methods
Synthesis of MoS 2 , WS 2 and WSe 2 . The MoS 2 monolayers were grown by low-pressure CVD technique for 10 min. The SiO 2 /Si substrates were cleaned using piranha solution (a volumetric mixture of 3:1 of 98% H 2 SO 4 to 35% H 2 O 2 ), then placed in the center of the quartz tube′ s heating zone and heated to 850 °C in argon atmosphere to restrict further oxidation. The molybdenum trioxide (MoO 3 ) powder was loaded in front of the substrates, while the sulfur powder was placed in the front of the quartz tube, which was heated to 200 °C by a heating band as shown in Figure S3. The WS 2 and WSe 2 monolayers were separately fabricated by constant-pressure vapor phase deposition method at 1100 °C and 1200 °C in argon atmosphere for 10 min.  Table 1. Fitting parameters of the PL peak energy as a function of temperature extracted in Fig. 4(f-h).
Scientific RepoRts | 7:44712 | DOI: 10.1038/srep44712 Transfer Method. The TMDs obtained by mechanical exfoliation have a small flake size, and it is not possible to create large scale heterostructured materials. Thus, the heterostructures here were prepared from chemical vapor deposition (CVD) growth (see Supporting Information Figure S3). Figure S4 illustrates the polystyrene (PS) film transfer technique schematically. The process started by coating a layer of polymer on top of the as-grown TMDs. The polymer acted as carrier layer, which can help with the handling of the atomically thin materials. The 10 wt% PS was dissolved in toluene. It was spin-coated on the CVD grown monolayer TMDs/SiO 2 /Si with a speed of 3500 rmp for 1 minute. Then, the samples with coating layer were baked at 90 °C for 1 hour. This was designed to evaporate the toluene and eliminate air bubbles formed at the interface of monolayer TMDs and PS, so that the adhesion is increased. This was followed by a gentle poking made by a glass cutter at the edge of the PS/TMDs/SiO 2 /Si to expose SiO 2 , by this way, the water penetration in the next step can be improved. Then, the PS/TMDs/SiO 2 /Si was sank into a large drop of deionized water in the hydrophobic plastic petri dishes gently. This water penetration operation was repeated for several times, until the PS/TMDs film was floated on the top of the water because of the repulsion between the PS and water. Next, the substrate was taken away, and another TMDs/SiO 2 /Si was put into the water from the side. After that, the system was dried by simply sucking the water droplet away with a paper towel. During this process, the PS/TMDs was fixed above the TMDs/SiO 2 /Si with a tweezers. After the PS/TMDs was transferred to the TMDs/SiO 2 /Si, the bubbles in the transferred assembly were purged out with the hydrogen. In order to combine the heterostructures better, the samples were baked at 90 °C for 30 min to remove water residues, and then at 120 °C for 30 min to spread the polymer for the elimination of possible wrinkles. Finally, the transferred assembly was soaked in toluene for at least 2 hours with the toluene changed several times during this time. It was worth noting that the obtained heterostructures were much cleaner than the conventional polymethyl methacrylate (PMMA) transfer techniques with reduced transfer film residue and improved quality ( Figure S5) 44,45 . Besides, the new transfer technique has the advantage of high efficiency on account of the regardless of the thickness of spin-coating layer.

Raman/PL Spectroscopy. Temperature dependent PL experiments were carried out by a Jobin-Yvon
LabRAM HR 800 micro-Raman spectrometer and a THMSE 600 heating/cooling stage (Linkam Scientific Instruments) in the temperature range from − 196 °C to 2 °C with a resolution of 0.1 °C. The heterostructures was excited by the 488 nm line of an Ar laser with the output power of 20 mW and recorded in back-scattering geometry with a resolution of better than 1 cm −1 . The laser beam was focused through a 50× microscope with a working distance of 18 mm. An air-cooled CCD (− 70 °C) with a 1024 × 256 pixels front illuminated chip was used to collect the scattered signal dispersed on 1800 grooves/mm grating for Raman and 600 grooves/mm grating for PL 46 . Peaks of the bare spectra are assigned by using the supporting software NGSLabSpec designed by JobinYvon.
Computational Details. Our density functional theory (DFT) calculations were performed with plane-wave pseudopotentials from the calculate method of quantum mechanics 47,48 . The generalized gradient approximation (GGA) for exchange-correlation term was employed, and the functional was adopted for structure optimization and band gaps calculation 49 . The selected pseudopotential is Ultrasoft Pseudopotentials. To acquire accurate results for bilayer structures, DFT-D approach was included with the Ortmann-Bechstedt-Schmidt (OBS) vdW correction 50 , which gives the structural parameters in good agreement with experimental values. A planewave cutoff energy of 320 eV and a 5 × 5 × 1 grid of Monkhorst-Pack points were employed to ensure good convergence of the computed structures and energies. Geometry optimization was determined using the Broyden-Fletcher-Goldfarb-Shenno (BFGS) minimization technique, with thresholds of converged structure of energy change per atom smaller than 10 −5 eV/atom, and displacement of atoms during geometry optimization no more than 0.001 A. The tolerance in the self-consistent field (SCF) calculation was 1.0 × 10 −6 eV/atom. All calculations were carried out using a 1 × 1 supercell with vacuum thickness not smaller than 17 Å and spin-orbit coupling was not included.