The article information
 Tianmin Wu, Ruixue Xu, Xiao Zheng, Wei Zhuang
 吴天敏, 徐瑞雪, 郑晓, 庄巍
 Electronic Structures and Thermoelectric Properties of TwoDimensional MoS_{2}/MoSe_{2} Heterostructures
 二维MoS_{2}/MoSe_{2}异质材料的电子结构和热电性质研究
 Chinese Journal of Chemical Physics , 2016, 29(4): 445452
 化学物理学报, 2016, 29(4): 445452
 http://dx.doi.org/10.1063/16740068/29/cjcp1512265

Article history
 Received on: December 30, 2015
 Accepted on: March 15, 2016
b. State Key Laboratory of Structural Chemistry, Fujian Institute of Research on the Structure of Matter, Chinese Academy of Sciences, Fuzhou 350002, China
Thermoelectric materials are considered to have great potential for power generation, energy saving, and environmental protection [15]. Various thermoelectric semiconductor materials, such as chalcogenides [6, 7], zintl phases [8], clatharates [9], complex oxides [10], and skutterudites [11, 12] have been developed to convert the waste heat into electricity. The efficiency of thermoelectric device is measured by the materialdependent figure of merit (ZT=S^{2}σT\κ) [13]. Here, S is the Seebeck coefficient, σ is the electrical conductivity, and κ=κ_{e}+κ_{lat} is the thermal conductivity (κ_{e}+κ_{lat} represent the electronic and lattice thermal conductivity contributions, respectively). In general, the figure of merit (ZT) can be enhanced by improving the power factor (S^{2}σ) or decreasing the thermal conductivity (κ) of material. The search of materials of high thermoelectric performance is confronted with the challenge that electronic conductivity and Seebeck coefficient have opposite dependencies on the materials parameters: electronic conductivity (σ) increases as the doping level is improved, but higher doping level leads to lower Seebeck coefficient. Various theoretical and experimental methods have been introduced to enhance the value of ZT [1426]. For instance, rattler atoms implanted into open structures like clatharates [27] and skutturides [28] have been introduced to have lowfrequency phonons near the acoustic branches, thus enhancing the phonon scattering to reduce the lattice thermal conductivity of the materials.
Twodimensional (2D) transitionmetal dichalcogenide semiconductors such as MoS_{2} and MoSe_{2} have significant potential as the ideal thermoelectric materials, since they possess low thermal conductivity along the direction perpendicular to the lattice layers [2933]. Despite this advantage, ZT values of transitionmetal dichalcogenides (TMDCs) are still low due to the difficulty of enhancing the power factor [3437]. To solve this problem, we propose the use of 2D heterostructures systems such as MoS_{2}/MoSe_{2}, which are formed by stacking two different TMDCs layers together. Based on density functional theory (DFT) calculations, we demonstrate that such 2D heterostructures materials possess many improved thermoelectric properties as compared with their parental pristine materials.
Ⅱ COMPUTATIONAL METHODSFirstprinciples calculations of bulk and bilayer MoS_{2}/MoSe_{2} 2D heterostructures system are conducted using DFT methods implemented in the Vienna ab initio simulation package (VASP) [38]. The projectoraugmentedwave pseudopotentials and the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) for exchangecorrelation functional are adopted in our simulations [39]. Furthermore, van der Waals interactions are taken into account by using the semiempirical correction of Grimme (DFTD3) [40]. The energy cutoff for a planewave expansion is set to 500 eV. All atomic coordinates are relaxed until the atomic forces have declined to 0.01 eV/Å, enforcing a total energy convergence criterion of 1×10^{5} eV. A vacuum slab larger than 15 Å is added to avoid interaction between adjacent images of the bilayer structure.
Temperatureand dopingdependent electronic transport properties, including electronic conductivity and Seebeck coefficient, are computed by using the semiclassical Boltzmann transport theory [41, 42]. The constant scattering time approximation is adopted, which is valid if the electron relaxation time does not very strongly with the energy on a scale of κ_{B}T, and the rigid band approaches as implemented in the BoltzTraP code [43]. By using a Fourier expansion [44, 45], while maintaining the crystal symmetry [46], the BoltzTrap code fits the ab initio electronic band structure to an analytic function. Since the transport properties can be very sensitive to the Brillouin zone (BZ) sampling, especially for low doping levels and low temperatures, we calculate the electronic structures required for the transport calculations with very dense kmeshes (43×43×11 for the hexagonal BZ of bulk 2D heterostructures, and 43×43×5 for the trigonal BZ of the bilayer structure).
The temperaturedependent and dopingleveldependent thermoelectric transport tensors, such as electronic conductivity σ_{αβ}(T, μ), Seebeck coefficient S_{αβ}(T, μ) and thermal conductivity (electronic part κ_{αβ}^{el}) tensors are calculated as follows [43]:
$ {\sigma _{\alpha \beta }}(T,\mu ) = \frac{1}{\Omega }\int {{{\bar \sigma }_{\alpha \beta }}} (\varepsilon )\left[{  \frac{{\partial {f_0}(T,\varepsilon ,\mu )}}{{\partial \varepsilon }}} \right]{\rm{d}}\varepsilon $  (1) 
$ \begin{array}{c} {S_{\alpha \beta }}(T,\mu ) = \frac{1}{{eT\Omega {\sigma _{\alpha \beta }}(T,\mu )}}\\ \int {{{\bar \sigma }_{\alpha \beta }}} (\varepsilon  \mu )\left[{  \frac{{\partial {f_0}(T,\varepsilon ,\mu )}}{{\partial \varepsilon }}} \right]{\rm{d}}\varepsilon \end{array} $  (2) 
$ \begin{array}{c} \kappa _{\alpha \beta }^{{\rm{el}}}(T,\mu ) = \frac{1}{{{e^2}T\Omega {\sigma _{\alpha \beta }}(T,\mu )}}\\ \int {{{\bar \sigma }_{\alpha \beta }}} {(\varepsilon  \mu )^2}\left[ {  \frac{{\partial {f_0}(T,\varepsilon ,\mu )}}{{\partial \varepsilon }}} \right]{\rm{d}}\varepsilon \end{array} $  (3) 
Here, α and β are the tensor indices, Ω is the volume of the unit cell, e is the charge of the electron, f_{0}(T, ε, μ) is the Fermi distribution function, and μ is the Fermi level. The conductivity tensor σαβ(ε) can be expressed analytically as [43]:
$ {\bar \sigma _{\alpha \beta }}(\varepsilon ) = \frac{1}{N}\sum\limits_{i,k} {{\sigma _{\alpha \beta }}} (i,{\bf{k}})\delta (\varepsilon  {\varepsilon _{i,k}}) $  (4) 
which can be expressed using kdepedent conductivity tensor as [43]:
$ {\sigma _{\alpha \beta }}(i,{\bf{k}}) = {e^2}{\tau _{i,k}}{\upsilon _\alpha }(i,{\bf{k}}){\upsilon _\beta }(i,{\bf{k}}) $  (5) 
where i is the band index, k is the reciprocal vector, N is a normalization depending on the number of k points sampled in the BZ, τ_{i, k} is the electronic relaxation time, and υα(i, k) is the i component of band velocity
The electron relaxation time τ characterizing the average time between two consecutive electron scattering events is a crucial parameter for calculating the thermoelectric properties, such as the electrical conductivity (σ) and the electronic thermal conductivity (κ_{e}), since just the ratio of these conductivity to the relaxation time (σ/τ and κ_{e}/τ) can be achieved by using the BoltzTrap code. The value of τ is usually obtained by fitting the calculated ratio of σ/τ to measure electrical conductivity data. However, such experimental data for the MoS_{2}/MoSe_{2} 2D heterostructures have so far remained unavailable. Therefore, in this work we focus on the ratio of thermal properties to the electron relaxation time rather than on the properties themselves.
Ⅲ RESULTS AND DISCUSSION A Crystal structureThe monolayer structure of MoS_{2} and MoSe_{2}, consists of a single (S/Se)Mo(S/Se) layer with space group P6m2 (187), which has no inversion symmetry. While the bulk MoS_{2} and MoSe_{2} has a 2Hpolytype structure in the P6_{3}/mmc space group (194). Its supercell consists of two (S/Se)Mo(S/Se) layers separated along the z axis, and the two layers are bound by van der Waals interactions. Furthermore, due to the increasing radius of the chalcogen atoms, the optimized values of lattice constants for monolayer along inplane direction increase from a=3.162 Å in MoS_{2} to 3.320 Å in MoSe_{2}. In general the lattice mismatch may lead to stacking disorder or Moiré Pattern superstructures. However, the intrinsic lattice mismatch between the layer of MoS_{2} and MoSe_{2} is as small as 0.158 Å. Moreover, the explicit consideration of such a small lattice mismatch would require the use of very large supercell and thus make the calculation rather expensive. Therefore, in our simulations for the 2D heterostructures, the same lattice constant is adopted for the both types of layers (MoS_{2} and MoSe_{2}).
The optimized MoS_{2}/MoSe_{2} 2D heterostructures for bulk system in our calculation has the lattice parameter of a=3.248 and 3.250 Å for bilayer, which is consistent with the other theoretical results [4749]. As in Fig. 1, the Mo atoms of MoSe_{2} monolayer sit on the top of the chalcogen atoms S of MoS_{2} monolayer, and those monolayer which constructs the heterostructures shows a lateral shift. Due to the van der Waals binding, the Mo atoms between different layers are separated by 6.36 Å in both bulk and bilayer systems. To further study the electronic properties and thermoelectric properties of bulk and bilayer MoS_{2}/MoSe_{2} heterostructures, their band structures and density of states (DOS) are calculated.
B Electronic structureElectronic band structures and DOS for bulk and bilayer MoS_{2}/MoSe_{2} 2D heterostructures are presented in Fig. 2. For the bulk 2D heterostructures, as demonstrated in Fig. 2(a), an indirect band gap of 0.546 eV is observed with the valence band maximum at the Γ point and the conduction band minimum at the K point, which is consistent with the theoretical results reported by Changhoon et al. [47]. The electronic band structure along high symmetry point of inplane and crossplane direction shows a great anisotropy due to the 2D heterostructures constructed by van der Waals binding, and the band gap of the inplane bands (ΓMKΓ) is much smaller than crossplane one (ΓA). Bilayer 2D heterostructures, on the other hand, also show an indirect band gap of 0.695 eV with both the valance band maximum at the Γ point and conduction band minimum at the high symmetry K point, which is consistent with other theoretical results [49]. Furthermore, both the bulk and bilayer systems show a strong asymmetric feature between the valence and conduction bands. Although GGAPBE is known to underestimate band gaps of semiconductors, the resulting electronic structure are considered to be reasonably accurate for subsequent computation of thermoelectric properties.
C Thermoelectric propertiesTo enlarge ZT, the material should have a larger power factor S^{2}σ and smaller κ (κ=κ_{e}+κ_{lat}. Figure 3 depict the inplane and crossplane power factor divided by the relaxation time τ, S^{2}σ/τ, of the bulk MoS_{2}/MoSe_{2} heterostructures under doping (both pand ntype), with the temperature ranging from 300 K to 1200 K. Here, τ is the relaxation time that is not directly determined by the band structure, but depends on the temperature, the doping level, and also the sample details (such as defect types and concentrations) [43]. Note that we compare S^{2}σ/τ instead of S^{2}σ, because the relaxation time is difficult to calculate and the electronic conductivity has not been measured experimentally. As demonstrated in Fig. 3, S^{2}σ/τ is enhanced as the carrier concentration increases at each temperature. Its maximum value is nearly within the carrier concentration range of 10^{20}10^{21} cm^{3} for each temperature, which also implies that high carrier concentration could enhance the power factor. As shown in Fig. 3, for bulk heterostructures, the ptype doping at inplane direction shows the largest power factor S^{2}σ/τ at temperature 1200 K. To better reflect the thermoelectric performance of bulk 2D MoS_{2}/MoSe_{2} heterostructures, the theoretical results of pristine bulk MoSe_{2} at the hole carrier concentration of 5×10^{20} cm^{3} are adopted [36]. As demonstrated in Table Ⅰ, although the pristine MoSe_{2} shows a slightly higher S^{2}σ/τ at 1200 K, the bulk MoS_{2}/MoSe_{2} heterostructures shows an overall better thermoelectric performance along the inplane direction. Despite the difference in relaxation time between these two crystals, such a comparison provides an overall assessment for their thermoelectric performances.
Along the crossplane direction, the ntype doping shows a greater power factor S^{2}σ/τ than ptype one. However, since in TMDCs the relaxation time of inplane direction is two orders of magnitude larger than the crossplane one [36, 37], the inplane electrical conductivity is typically two orders magnitude larger than the crossplane electrical conductivity. Then, the power factor along inplane direction is expected to be two orders of magnitude larger than the crossplane counterpart, since both of them have similar Seebeck coefficient values. Additionally, the DOS close to the valence band edge is much larger than those near the conduction band edge, which suggests that ptype doping could have a better thermoelectric performance [50]. Therefore, we will focus on the ptype doping bulk heterostructures, and discuss the temperaturedependent and dopingleveldependent Seebeck coefficient, the electric conductivity and the thermal conductivity individually.
The ptype doping inplane and crossplane Seebeck coefficients as functions of carrier concentration are shown in Fig. 4 (a) and (b). Similar to the other TMDCs [36, 37], the bulk 2D heterostructures has a large Seebeck coefficient. At a fixed carrier concentration, for both inplane and crossplane directions, the Seebeck coefficient of bulk 2D heterostructures slightly increases as the temperature is increasing. This phenomenon mainly originates from the Fermi broadening as the temperature rises, and then it leads to an increasing effective density of states at the top of valence band [51]. Consistent with the known thermoelectric behavior for the other TMDCs [36, 37], the maximum value of Seebeck coefficient for each temperature shifts to high doping level and decrease as temperature is increased for both crossand inplane direction. The calculated tiny band gap for bulk 2D heterostructures (0.546 eV), as mentioned above, likely lead to the bipolar effect at low doping carrier concentration which make the Seebeck coefficient decreases with decreasing doping concentration, opposite to the usual situation [52]. As temperature is higher than 900 K, the Seebeck coefficient along the inplane direction presents a sign reversal at low doping level, which is attributed to the increasing negative contribution of thermally excited electrons to the Seebeck coefficient under the bipolartransport conditions [53]. While the Seebeck coefficient along the crossplane direction does not show a sign reversal and more stable with the doping concentration increasing at low doping level, indicating that the system along crossplane direction shows a weak bipolar effect originated by the large band gap along this direction. This anisotropy can also be observed by the electrical conductivity divided by the relaxation time (σ/τ) as a function of doping concentration. As demonstrated in Fig. 4 (c) and (d), the conductivity divided by the relaxation time (σ/τ) of inplane is about one order higher than the crossplane one, revealing that the thermally excited carrier concentration along the inplane direction is significantly higher than crossplane one. Furthermore, in the 9001200 K temperature range, the electrical conductivity divided by the relaxation time (σ/τ) does not depend strongly on the doping concentration at low doping level, since the thermally excited carriers dominate transport [53]. To overcome the contribution by the thermally excited carrier concentration, the higher doping concentration is required as temperature is higher. As the doping concentration above the onset of bipolar transport, the electrical conductivity divided by the relaxation time (σ/τ) increases substantially with the increasing doping concentration.
The electronic band structure also supports this viewpoint, as the band gap of inplane (along high symmetry point ΓMKΓ) is much smaller than the crossplane one (ΓA). The electrical thermal conductivities divided by the relaxation time (κ_{e}/τ) along the inplane and crossplane directions are depicted in Fig. 4 (e) and (f), respectively. Although the electronic contribution to the thermal conductivity is small, compared with the lattice thermal conductivity, the electronic thermal conductivity also illustrates the anisotropy of thermal conductivity between crossplane and inplane. The power factor of inplane is considerably larger than the crossplane one, especially since the relaxation time of crossplane is much smaller than the inplane one. Furthermore, the electronic thermal conductivity divided by the relaxation time (κ_{e}/τ) at a fixed carrier concentrations for both inplane and crossplane increases as temperature rising, mainly because of more intense electron scattering as temperature increasing. According to the WiedemannFranz law [43],
$ \kappa _{\alpha \beta }^{{\rm{el}}} = \frac{{{\pi ^2}}}{3}{\left( {\frac{{{k_{\rm{B}}}}}{e}} \right)^2}{\sigma _{\alpha \beta }}T $  (6) 
the electronic thermal conductivity is proportional to the temperature of system. Meanwhile, since thermal conductivity is proportional to the electronic conductivity that is highly anisotropic, there is a significant difference between the inplane and crossplane components.
Compared with the bulk heterostructures, the bilayer heterostructures show a much lower power factor (S^{2}σ/τ) as demonstrated in Fig. 5. Since the band gap of bilayer 2D heterostructures is larger than the bulk one, the electronic conductivity of bilayer heterostructures is much smaller than that of the bulk, see Fig. 6 (a) and (b). Furthermore, since the bulk and bilayer heterostructures have similar Seebeck coefficients, as demonstrated in Fig. 6 (c) and (d), the bulk heterostructures should have a better thermoelectric performance than those of bilayer one. While the electrical thermal conductivity is smaller than the bulk one, as presented in Fig. 6 (e) and (f). For ptype doping, the bilayer heterostructures show a similar thermoelectric performance to the monolayer one [36].
Ⅳ CONCLUSIONBased on the electrical band structure calculated from first principles, the thermoelectric properties of bulk and bilayer 2D MoS_{2}/MoSe_{2} heterostructures have been analyzed by using the semiclassical Boltzmann transport theory. Both ntype and ptype doping have been addressed for bulk and bilayer heterostructures, employing the rigid band approximation and constant scattering time approximation. Due to a smaller band gap and more dense DOS close to the valence band edge than those near the conduction band edge, the thermoelectric performance of bulk 2D MoS_{2}/MoSe_{2} heterostructures turns out to be superior to pristine bulk MoSe_{2} along the inplane direction for ptype doping at a wide temperature range. Furthermore, with a larger band gap, the bilayer heterostructures show a much lower electronic conductivity than those of bulk heterostructures, which also induces that it shows a weaker thermoelectric performance than bulk one. Although the power factor of bilayer heterostructures is lower than those of bulk heterostructures, it shows a similar thermoelectric performance to the monolayer MoSe_{2} for ptype doping at each temperature. Therefore, we safely conclude that such 2D heterostructures materials possess much improved thermoelectric properties as compared with their parental pristine materials, especially for bulk one. As reported by Li et al. [54], the thermal lattice conductivity considered the spinorbit coupling (SOC) effect is higher than those value not considered. Furthermore, the band shape and the band gap which make a great influence on the calculation of thermoelectric properties could be changed as the SOC effect is introduced, thus the SOC effect should be taken into account. However, the SOC is not considered in the present work, because of the limited computational resources at our disposal. The influence of the SOC on the thermoelectric properties is to be addressed in our future work.
Ⅴ ACKNOWLEDGMENTSThis work was supported by the National Natural Science Foundation of China (No.21203178, No.21373201, No.21433014, No.21233007, No.21303175, and No.21322305), the Science and Technological Ministry of China (No.2011YQ09000505), the ''Strategic Priority Research Program'' of the Chinese Academy of Sciences (No.XDB10040304 and No.XDB100202002), and the Fundamental Research Funds for the Central Universities (No.2340000074). The computational resources are provided by the Supercomputing Center of University of Science and Technology of China.
[1]  M. S. Dresselhaus, G. Chen, M. Y. Tang, R. G. Yang, H. Lee, D. Z. Wang, Z. F. Ren, J. P. Fleurial,and P. Gogna, Adv. Mater. 19 , 1043 (2007). DOI:10.1002/(ISSN)15214095 
[2]  K. Biswas, J. He, I. D. Blum, C. Wu, T. Hogan, D. N. Seidman, V. P. Dravid,and M. G. Kanatzidis, Nature 489 , 414 (2012). DOI:10.1038/nature11439 
[3]  K. P. Pernstich, B. Rössner,and B. Batlogg, Nat. Mater. 7 , 321 (2008). DOI:10.1038/nmat2120 
[4]  W. J. Liang, A. I. Hochbaum, M. Fardy, O. Rabin, M. J. Zhang,and P. D. Yang, Nano Lett. 9 , 1689 (2009). DOI:10.1021/nl900377e 
[5]  G. J. Snyder,and E. S. Toberer, Nat. Mater. 7 , 105 (2008). DOI:10.1038/nmat2090 
[6]  C. Wood, Rep. Prog. Phys. 51 , 459 (1988). DOI:10.1088/00344885/51/4/001 
[7]  L. E. Shelimova, O. G. Karpinskii, P. P. Konstantinov, E. S. Avilov, M. A. Kretova,and V. S. Zemskov, Inorg. Mater. 40 , 451 (2004). DOI:10.1023/B:INMA.0000027590.43038.a8 
[8]  S. M. Kauzlarich, S. R. Brown,and G. J. Snyder, Dalton Trans. , 2009 (2007). 
[9]  G. S. Nolas, G. A. Slack, and S. B. Schujman, In Recent Trends in Thermoelectric Materials Research I, T. M. Tritt, Ed., Semiconductors and Semimetals, San Diego, CA: Academic Press, 255 (2001). 
[10]  K. Koumoto, I. Terasaki,and R. Funahashi, MRS Bull. 31 , 206 (2006). DOI:10.1557/mrs2006.46 
[11]  G. S. Nolas, D. T. Morelli,and T. M. Tritt, Annu. Rev. Mater. Sci. 29 , 89 (1999). DOI:10.1146/annurev.matsci.29.1.89 
[12]  C. Uher, In Recent Trends in Thermoelectric Materials Research I, T. M. Tritt, Ed., Semiconductors and Semimetals, San Diego, CA: Academic Press, 139 (2001). 
[13]  H. J. Goldsmid, In Recent Trends in Thermoelectric Materials Research I, T. M. Tritt, Ed., Semiconductors and Semimetals, San Diego, CA: Academic Press, 1 (2001). 
[14]  G. A. Slack,and D. M. Rowe, CRC Handbook of Thermoelectrics. Boca Raton: CRC Press, 40 (1995). , 40 (1995). 
[15]  L. D. Hicks,and M. S. Dresselhaus, Phys. Rev. B 47 , 12727 (1993). DOI:10.1103/PhysRevB.47.12727 
[16]  L. D. Hicks,and M. S. Dresselhaus, Phys. Rev. B 47 , 16631 (1993). DOI:10.1103/PhysRevB.47.16631 
[17]  Y. Wu, R. Fan,and P. Yang, Nano Lett. 2 , 83 (2002). DOI:10.1021/nl0156888 
[18]  R. Yang, G. Chen,and M. S. Dresselhaus, Phys. Rev. B 72 , 125418 (2005). DOI:10.1103/PhysRevB.72.125418 
[19]  C. J. Vineis, A. Shakouri, A. Majumdar,and M. G. Kanatzidis, Adv. Mater. 22 , 3970 (2010). DOI:10.1002/adma.201000839 
[20]  A. I. Boukai, Y. Bunimovich, J. TahirKheli, J. K. Yu, W. A. Goddard Ⅲ,and J. R. Heath, Nature 451 , 168 (2008). DOI:10.1038/nature06458 
[21]  X. Tang, W. Xie, H. Li, W. Zhao, Q. Zhang,and M. Niino, Appl. Phys. Lett. 90 , 012102 (2007). DOI:10.1063/1.2425007 
[22]  B. Poudel, Q. Hao, Y. Ma, Y. Lan, A. Minnich, B. Yu, X. Yan, D. Wang, A. Muto, D. Vashaee, X. Chen, J. Liu, M. S. Dresselhaus, G. Chen,and Z. Ren, Science 320 , 634 (2008). DOI:10.1126/science.1156446 
[23]  T. Ikeda, L. A. Collins, V. A. Ravi, F. S. Gascoin, S. M. Haile,and G. J. Snyder, Chem. Mater. 19 , 763 (2007). DOI:10.1021/cm062121p 
[24]  R. Venkatasubramanian, E. Siivola,and T. Colpitts, B. O'Quinn, Nature 413 , 597 (2001). 
[25]  T. Ikeda, S. M. Haile, V. A. Ravi, H. Azizgolshani, F. Gascoin,and G. J. Snyder, Acta Mater. 55 , 1227 (2007). DOI:10.1016/j.actamat.2006.09.036 
[26]  T. C. Harman, P. J. Taylor, M. P. Walsh,and B. E. LaForge, Science 297 , 2229 (2002). DOI:10.1126/science.1072886 
[27]  J. Dong, O. F. Sankey,and C. W. Myles, Phys. Rev. Lett. 86 , 2361 (2001). DOI:10.1103/PhysRevLett.86.2361 
[28]  V. Keppens, D. Mandrus, B. C. Sales, B. C. Chakoumakos, P. Dai, R. Coldea, M. B. Maple, D. A. Gajewski, E. J. Freeman,and S. Bennington, Nature 395 , 876 (1998). DOI:10.1038/27625 
[29]  L. H. Brixner, J. Inorg. Nucl. Chem. 24 , 257 (1962). DOI:10.1016/00221902(62)80178X 
[30]  S. H. El,  Mahalawy,and B. L. Evans, Phys. Status Solidi B 79 , 713 (1977). DOI:10.1002/(ISSN)15213951 
[31]  E. Revolinsky,and D. Beerntsen, J. Appl. Phys. 35 , 2086 (1964). DOI:10.1063/1.1702795 
[32]  C. Muratore, V. Varshney, J. J. Gengler, J. J. Hu, J. E. Bultman, T. M. Smith, P. J. Shamberger, B. Qiu, X. Ruan, A. K. Roy,and A. A. Voevodin, Appl. Phys. Lett. 102 , 081604 (2013). DOI:10.1063/1.4793203 
[33]  A. Mavrokefalos, N. T. Nguyen, M. T. Pettes, D. C. Johnson,and L. Shi, Appl. Phys. Lett. 91 , 171912 (2007). DOI:10.1063/1.2800888 
[34]  M. Kayyalha, L. Shi, and Y. P. Chen, arXiv: 1505.05891 (2015). 
[35]  H. Kedar, W. Ying, Y. Yu, Z. Hanyu, W. Yuan, M. Joel, and Z. Xiang, arXiv:1505.06779 (2015). 
[36]  S. Kumar,and U. Schwingenschlöogl, Chem. Mater. 27 , 1278 (2015). DOI:10.1021/cm504244b 
[37]  A. N. Gandi,and U. Schwingenschlöogl, Chem. Mater. 26 , 6628 (2014). DOI:10.1021/cm503487n 
[38]  G. Kresse,and J. Furthmüller, Phys. Rev. B 54 , 11169 (1996). DOI:10.1103/PhysRevB.54.11169 
[39]  J. P. Perdew, K. Burke,and M. Ernzerhof, Phys. Rev. Lett. 77 , 3865 (1996). DOI:10.1103/PhysRevLett.77.3865 
[40]  S. Grimme,and G. G. Semiempirical, J. Comput. Chem. 27 , 1787 (2006). DOI:10.1002/(ISSN)1096987X 
[41]  J. M. Ziman, Principles of the Theory of Solids (1972). 
[42]  P. B. Allen, W. E. Pickett,and H. Krakauer, Phys. Rev. B 37 , 7482 (1988). DOI:10.1103/PhysRevB.37.7482 
[43]  G. Madsen,and D. Singh, Comput. Phys. Commun. 175 , 67 (2006). DOI:10.1016/j.cpc.2006.03.007 
[44]  R. N. Euwema, D. J. Stukel, T. C. Collins, J. S. Dewitt,and D. G. Shankland, Phys. Rev. 178 , 1419 (1969). DOI:10.1103/PhysRev.178.1419 
[45]  W. Jones,and N. March, Theoretical Solid State Physics: Perfect Lattices in Equilibrium, Dover Books on Physics (1973). 
[46]  D. D. Koelling,and J. H. Wood, J. Comput. Phys. 67 , 253 (1986). DOI:10.1016/00219991(86)902615 
[47]  L. Changhoon, H. Jisook, W. MyungHwan,and H. S. Ji, Chem. Mater. 25 , 3745 (2013). DOI:10.1021/cm402281n 
[48]  T. Humberto, L. Florentino,and T. Mauricio, Sci. Rep. 3 , 1549 (2013). 
[49]  N. Lu, H. Guo, L. Lei, J. Dai, L. Wang, W. N. Mei, X. Wu,and X. C. Zeng, Nanoscale 6 , 2879 (2014). DOI:10.1039/c3nr06072a 
[50]  Z. Lijun,and J. S. David, Phys. Rev. B 81 , 245119 (2010). DOI:10.1103/PhysRevB.81.245119 
[51]  Z. Gang,and W. Dong, Sci. Rep. 5 , 8099 (2015). DOI:10.1038/srep08099 
[52]  S. BaoZhen, M. Zuju, H. Chao,and W. Kechen, RSC Adv. 5 , 56382 (2015). DOI:10.1039/C5RA06134B 
[53]  G. Shi,and E. Kioupakis, J. Appl. Phys. 117 , 065103 (2015). DOI:10.1063/1.4907805 
[54]  R. Li, X. Cheng, Q. Xie, Y. Sun, D. Li, Y. Li,and X. Chen, Sci. Rep. 5 , 8466 (2015). DOI:10.1038/srep08466 
b. 中国科学院福建物质结构研究所, 结构化学国家重点实验室, 福州 350002