Bonding Charge Density and Ultimate Strength of Monolayer Transition Metal Dichalcogenides
BBonding Charge Density and Ultimate Strength of Monolayer Transition MetalDichalcogenides
Junwen Li, Nikhil V. Medhekar, and Vivek B. Shenoy ∗ Department of Materials Science and Engineering,University of Pennsylvania, Philadelphia, PA 19104, USA Department of Materials Engineering, Monash University, Clayton, Victoria 3800, Australia (Dated: April 1, 2013)Two-dimensional (2D) semiconducting transition metal dichalcogenides (TMDs) can withstand alarge deformation without fracture or inelastic relaxation, making them attractive for applicationin novel strain-engineered and flexible electronic and optoelectronic devices. In this study, wecharacterize the mechanical response of monolayer group VI TMDs to large elastic deformationusing first-principles density functional theory calculations. We find that the ultimate strength andthe overall stress response of these 2D materials is strongly influenced by their chemical compositionand loading direction. We demonstrate that differences in the observed mechanical behavior canbe attributed to the spatial redistribution of the occupied hybridized electronic states in the regionbetween the transition metal atom and the chalcogens. In spite of the strong covalent bondingbetween the transition metal and the chalcogens, we find that a simple linear relationship can beestablished to describe the dependence of the mechanical strength on the charge transfer from thetransition metal atom to the chalcogens.
I. INTRODUCTION
Since its first mechanical exfoliation from graphite,graphene — an allotrope of carbon with a two-dimensional (2D) structure — has attracted tremendousscientific and technological attention due to novel elec-tronic, mechanical, and chemical properties.
The in-tense interest in graphene has also stimulated an activesearch for possible inorganic 2D or quasi-2D materialswith unique characteristics. Examples of such mate-rials include hexagonal boron nitride, transition metaldichalcogenides (TMDs), and transition metal oxides.
Earlier attempts to obtain monolayer and a few-layerthick samples of these materials were based on mechani-cal exfoliation from their bulk counterparts. Recent ad-vances in fabrication techniques, such as liquid exfolia-tion, lithium intercalation, epitaxy, laser thinning, andchemical vapor deposition, have successfully resulted ina scalable production of large 2D samples containing onlya single layer to a few layers.
Among the non-carbon 2D materials that are beingexplored, dichalcogenides of group IV, V and VI transi-tion metals are of particular interest. Similar to grapheneand its bulk counterpart graphite, these materials havestrong in-plane bonding while the individual layers arebonded by weak van der Waals interactions. SeveralTMDs, unlike graphene, have an intrinsic band gap in therange of 1 - 2 eV, making them attractive for field effecttransistors and optoelectronic devices.
Many single-layer TMDs also demonstrate novel electronic propertiesthat are otherwise not seen in their bulk or few-layerthick counterparts, for example, direct to indirect bandgap transition, enhanced photoluminescence, and valleypolarization.
Recent mechanical experiments show that 2D TMDscan sustain very large elastic strains ( ∼
10% effective in-plane strain) with a high resistance to inelastic relax-ation and fracture, and hence these materials are also ofinterest for flexible electronics applications. Since theelectronic conduction and valence energy states both de-pend on the strain, it has also been proposed that elasticstrain can be employed to tune the band gap in mono-layer TMDs.
In order to realize the potential ap-plications of various 2D TMDs in strain-engineered andflexible devices, an accurate characterization of their abil-ity to withstand large elastic deformation is crucial.Recently, Bertolazzi et al. performed the experimen-tal measurements on the stiffness and breaking strengthof monolayer MoS by using an atomic force micro-scope tip to deform the monolayer MoS placed on aprepatterned SiO substrate. They found that the ef-fective Young’s modulus and average breaking strengthof monolayer MoS are 270 ±
100 GPa and 23 GPa, re-spectively. Using a similar technique, Castellanos-Gomez et al. studied the mechanical properties of freely sus-pended MoS nanosheets with 5 to 25 layers and ob-tained the mean Young’s modulus of 330 ±
70 GPa.These experiments clearly establish 2D MoS as an ultra-strength material able to withstand large elastic defor-mation. However, such studies provide only effectiveorientation-averaged values of the mechanical propertiesand are limited only to single composition, while the restof the 2D semiconducting TMDs remain yet to be charac-terized. The theoretical studies of mechanical propertiesof monolayer TMDs are also limited to the calculation ofsimple elastic properties that are relevant only at smalldeformations. What is the intrinsic ultimate strength— the maximum stress a defect-free material can with-stand without failure — of monolayer TMDs along dif-ferent crystal directions and how is it dependent on thechemical composition? Indeed, the relationship betweenthe intrinsic stress response of monolayer TMDs to largeelastic deformation, the atomic-level geometry, and the a r X i v : . [ c ond - m a t . m t r l - s c i ] M a r chemical composition remains unclear.In this article, using first-principles density functionaltheory simulations, we report on the fundamental me-chanical behavior of monolayer TMDs when subjected toa large deformation. Here we consider TMDs with chem-ical composition MX , where M (= Mo, W) is a groupVI transition metal and X (= S, Se, Te) is a chalcogen.Among the range of naturally occurring TMDs, group VITMDS are semiconductors and, therefore, are good can-didates for flexible and strain-engineered optoelectronicand photonic devices. Here we have determined the en-tire stress response of monolayer TMDs to large appliedelastic strains — the intrinsic materials characteristicsthat are crucial in designing the strain-engineered andflexible devices. We find that the stress response andthe ultimate strength of monolayer TMDs strongly de-pends on the chemical composition as well as the load-ing direction. For the same transition metal M, the sul-fides (MS ) are the strongest and selenides (MSe ) arestronger than the tellurides (MTe ). The chalcogens ofW (WX ) can accommodate larger stresses than that ofMo (MoX ). For all compositions, the ultimate strengthof the monolayer along the armchair direction far exceedsthe strength in zigzag direction, typically by a factor of1.5 - 2.2. Our electronic structure analysis reveals thatthe origin of the observed mechanical behavior can beattributed to the hybridization between M- d and X- p or-bitals and the charge transfer from the transition metalto the chalcogens. II. RESULTS AND DISCUSSION
All MX TMDs monolayers considered in our studyshare a similar hexagonal crystal structure as shown inFig. 1. A single layer of MX consists of three X-M-X sublayers, with hexagonal lattice plane of the transi-tion metal M sandwiched between identical hexagonalplanes of the chalcogens X. In order to evaluate themechanical response of the pristine TMD monolayers,we carried out first-principles density functional calcu-lations within the local density approximation as im-plemented in ABINIT. Norm-conserving pseudopo-tentials in the form of Hartwigsen-Geodecker-Hutterpseudopotentials including the semicore electrons asvalence electrons for Mo and W were used to describethe interaction between core and valence electrons. Anenergy cutoff of 65 Hartree was used for plane wave basisexpansion. The k -point sampling schemes of 15 × × × × monolayers, we first performedthe structure relaxation including the lattice constantsand atomic coordinates as summarized in Table I. Thelattice constants obtained from our DFT calculations arein excellent agreement with the experimental measure-ments (error less than 1%). Our results are also in close x (a)(b) y FIG. 1. A schematic showing the crystal structure of a TMDmonolayer (MX ) with (a) top view and (b) side view. Thetransition metal (M) and chalcogen (X) atoms are representedby purple (large) and yellow (small) spheres, respectively. Theunit cell (red) and the orthogonal supercell (blue) are alsodepicted in (a). The armchair and zigzag orientations corre-spond to x and y directions, respectively.TABLE I. Lattice constants, M-X bond lengths and X-X dis-tances for monolayer MX (M = Mo, W; X = Se, Se, Te)TMDs.MX a (˚A) M-X (˚A) X-X (˚A)MoS agreement with those computed with different methodssuch as numerical atomic orbitals, PAW potentials .Similar to the 2D hexagonal lattice of graphene, twohigh-symmetry directions can be identified in the crys-tal structure of monolayer TMDs, namely, armchair andzigzag directions, which are oriented along and perpen-dicular to M-X bonds when projected on the plane con-taining transition metal atoms, respectively. Since themechanical behavior of a defect-free material is ulti-mately controlled by the strength of its chemical bonds, itcan be expected that the stress response of the monolay-ers along the armchair and zigzag directions will also bedistinct. We therefore employed an orthogonal supercellas shown in Fig. 1(a) to calculate the stress-strain curvesalong the two high symmetry directions. In each case, thesuper-cell was deformed along the high symmetry direc-tion in small strain increments and the in-plane supercellvector normal to the the applied strain was allowed to xy xy (b)(a) FIG. 2. Tensile stress σ as a function of uniaxial strain ε along (a) armchair and (b) zigzag directions, respectively, for monolayerMX (M = Mo, W; X = S, Se, Te) TMDs. Solid and dashed lines are used for WX and MoX , respectively. relax in order to account for the Poisson’s contractionas typically observed in the experiments. The supercellstress as output directly by ABINIT is the stress aver-aged over the supercell volume. We renormalized thecalculated supercell stress by a factor of Z/h , where Z isthe vacuum distance 12 ˚A used to avoid the interactionbetween periodic images and h is the interlayer distancein the bulk MX . The strain is defined as ε = L − L L ,where L and L are the supercell length along armchair orzigzag directions before and after applying tensile strain,respectively.In Figs. 2(a) and 2(b) we present the calculated stress-strain relations along armchair ( x ) and zigzag ( y ) direc-tions, respectively. It can be readily noted that for smallstrain, the stress for all MX exhibits linear dependenceon the applied strain for both loading directions. Ta-ble II presents the Young’s modulus E calculated as theslope of the stress-strain curve in the region of strain lessthan 4%. It is evident that for all compositions, the val-ues of the Young’s modulus obtained independently fromthe stress-strain curve for armchair and zigzag loadingdirection are virtually identical. This behavior is consis-tent with the symmetry of the MX crystal lattice - theYoung’s modulus and other second order elastic constantsare essentially isotropic due to the hexagonal symmetryin the basal plane. Our calculated values for the Young’smodulus are also in good agreement with earlier studies.For instance, we obtained the Young’s modulus of MoS to be ∼
220 GPa, while experimental measurements havereported values of 270 ±
100 GPa and 330 ±
70 GPa formonolayer and few layer MoS sheets , respectively.As MX monolayers are strained further ( ε > σ ∗ . Table II presents the cal-culated values for the peak stress σ ∗ and the corre-sponding ultimate strain ε ∗ for all chemical composi-tions. It can be observed that in general, the chalcogensof W (WX ) have larger moduli and tensile strength thanthat of Mo (MoX ). For the same transition metal, sul-phides (MS ) are the strongest while tellurides (MTe )are weakest. Moreover, Table II also presents valuesfor the anisotropy factor φ = σ ∗ AR /σ ∗ ZZ , which measuresthe relative strength along the two high symmetry direc-tions. We find that the anisotropy in stress response isinversely correlated with the strength of the monolayersheets — the MX with lower Young’s modulus and ulti-mate strength (for example, tellurides) are characterizedby larger anisotropy factors.Our calculated value for the average ultimate strengthof monolayer MoS (22 GPa) is in excellent agreementwith the value of 23 GPa obtained by the atomic forcemicroscopy measurements. Furthermore, the observedstress-strain response of the single-layer structure ofthe TMDs can be compared with their nanotube coun-terparts. For instance, density functional tight bind-ing (DFTB) simulations have predicted a Young’s modu-lus of 209.7 GPa and 236.6 GPa for armchair and zigzagnanotubes of comparable diameter about 12.0 ˚A, respec-tively, in close agreement with our results. Moreover,their ultimate strength was reported to be 29.1 GPa and31.6 GPa, respectively, again in good qualitative agree-ment with our results. However the DFTB studies ofMoS nanotubes suggested similar value ( ∼ ε ∗ corresponding to the ultimate strengthfor both armchair and zigzag directions, in contrast toour results where much bigger values of ε ∗ ( ≥ ε ∗ ( ≤ TABLE II. Calculated Young’s modulus E (GPa), ultimate strength σ ∗ (GPa), ultimate strain ε ∗ corresponding to the ultimatestrength for armchair and zigzag directions, anisotropy factor φ and the charge transfer ∆ Q ( e ) from M to X for strain-freeMX monolayer sheets.MX Armchair ( x ) Zigzag ( y ) φ ∆ QE σ ∗ ε ∗ E σ ∗ ε ∗ MoS MoS TDOS . . . . . D e n s i t y o f S t a t e s ( S t a t e s / e V U n i t C e ll ) Mo- s Mo- p Mo- d − − − Energy (eV) . . . . . . S- s S- p (a) MoSe TDOS . . . . . D e n s i t y o f S t a t e s ( S t a t e s / e V U n i t C e ll ) Mo- s Mo- p Mo- d − − − Energy (eV) . . . . . . Se- s Se- p (b) MoTe TDOS . . . . . D e n s i t y o f S t a t e s ( S t a t e s / e V U n i t C e ll ) Mo- s Mo- p Mo- d − − − Energy (eV) . . . . . . Te- s Te- p (c) WS TDOS . . . . . D e n s i t y o f S t a t e s ( S t a t e s / e V U n i t C e ll ) W- s W- p W- d − − − − Energy (eV) . . . . . . S- s S- p (d) WSe TDOS . . . . . D e n s i t y o f S t a t e s ( S t a t e s / e V U n i t C e ll ) W- s W- p W- d − − − − Energy (eV) . . . . . . Se- s Se- p (e) WTe TDOS . . . . . D e n s i t y o f S t a t e s ( S t a t e s / e V U n i t C e ll ) W- s W- p W- d − − − Energy (eV) . . . . . . Te- s Te- p (f) FIG. 3. Total and projected density of states for strain-free monolayer MX (M = Mo, W; X = S, Se, Te) TMDs. For thetransition metal M and the chalcogen X, the density of states are projected onto s , p , d and s , p orbitals, respectively. TheFermi level is set to zero. response at large strains has also been observed in thecase of graphene. However, there is one remarkable dif-ference: while TMDs and graphene both have hexagonallattice, they exhibit a contrasting anisotropy in stress re-sponse. The peak stress in graphene for armchair tensionis 9% lower than that for zigzag direction, in contrast tothe behavior seen in Fig. 2.In order to identify the atomic-level origin of the me-chanical behavior of various transition metal dichalco- genides considered here, we analyzed their electronicstructures by using the Vienna ab initio simulation pack-age (VASP) with LDA exchange correlation functional.Frozen-core projector-augmented wave (PAW) methodwas employed to describe the interaction between coreelectrons and valence electrons. An energy cutoff of 400eV was used for plane wave basis expansion.The calculated total and projected density of states ofunstrained MX monolayer sheets are presented in Fig. 3.Since the strength of the interatomic bonding and theconsequent mechanical response of the materials to defor-mation are determined primarily by the occupied statesjust below Fermi level, we focus on the energy range 0 -8 eV below the Fermi level. It can be seen that in thisrange, the outermost d orbitals of the transition metal Moverlap significantly with the outermost p orbitals of thechalcogen X, indicating a strong p - d hybridization. Fur-thermore, a careful examination of the peaks of these hy-bridized states yields a crucial insight on the compositiondependence of the mechanical response seen in Fig. 2. Forinstance, the peaks of p - d orbitals in MoS are locatedat -5.25, -4.02, -3.27, -1.93 and -0.83 eV below the Fermilevel. In MoSe , these peaks shift toward the Fermi levelby 0.28, 0.13, 0.15, 0.03 and 0.11 eV, respectively. Thepeaks shift further toward the Fermi level in MoTe by0.31, 0.52, 0.78, 0.30 and 0.08 eV, respectively. In addi-tion to the change in peak positions, the height of peaksalso changes as the composition is changed from MoS to MoTe : the peaks of low-lying Mo-4 d states are low-ered while the peaks of states close to Fermi level areenhanced. As seen in Fig. 3, the hybridized states forthe chalcogens of W follow an essentially similar trend.The changes in the p - d hybridized orbitals in monolayerTMDs can be quantitatively characterized by computingthe center of d bands below Fermi level, which is definedas E d = (cid:82) E F E L PDOS(
E, d ) × E dE (cid:82) E F E L PDOS(
E, d ) dE where PDOS( E, d ) is the density of states projected ontothe d orbitals of transition metals Mo and W, the lowerend energy E L is taken as -8 eV to cover the p - d hy-bridized states and the Fermi level E F is set to zero.For MoS , MoSe and MoTe , the calculated d bandcenters E d are -2.68, -2.37, -1.94 eV, respectively. ForWS , WSe and WTe , E d is located at -2.98, -2.65, -2.18 eV, respectively. Since the hybridization of p - d or-bitals indicate a covalent electron sharing between M andX atoms, the deeper the center of the d bands below theFermi level, the stronger is the M-X bond. Consequently,the tellurides with shallow E d demonstrate lower Young’smodulus and tensile strength compared to sulphides andselenides. Chalcogens of W have deeper E d than thechalcogens of Mo and are therefore stronger. The analysisof the locations of d band centers presented here explainsthe observed composition dependence of the mechanicalresponse of the monolayer TMDs as seen in Fig. 2.In order to visualize the hybridized electronic states be-tween the transition metals and chalcogens, we have fur-ther calculated the electronic charge distribution. Fig. 4shows the bonding charge density in the plane passingthrough both the transition metal and chalcogen atoms.The bonding charge density is obtained as the differencebetween the valence charge density of strain-free MX sheet and the superposition of the valence charge densityof the constituent atoms. A positive value (red) indicates WS MoS MoSe MoTe WSe WTe (a) (c) (e)(b) (d) (f) -0.01 FIG. 4. Bonding charge density for monolayer MX (M =Mo, W; X = S, Se, Te) TMDs, obtained as the charge densitydifference between the valence charge density of the monolayerand the superposition of the valence charge density of theneutral constituent atoms. Red and blue colors indicate theelectron accumulation and depletion, respectively. The colorscale is in the units of e /Bohr . . . . . . Charge transfer ∆ Q ( e ) Y o un g ’ s m o d u l u s E ( G P a ) WS MoS WSe MoSe WTe MoTe ArmchairZigzag (a) . . . . . Charge transfer ∆ Q ( e ) U l t i m a t e s t r e n g t h σ ∗ ( G P a ) WS MoS WSe MoSe WTe MoTe ArmchairZigzag (b)
FIG. 5. Variation of (a) the Young’s modulus and (b) theultimate strength of monolayer MX (M = Mo, W; X = S,Se, Te) TMDs with the charge transfer ∆ Q from transitionmetal M to chalcogens X. TABLE III. Fitted values in the linear relation ∼ a + b ∆ Q to describe the dependence of the Young’s modulus E andultimate strength σ ∗ AR and σ ∗ ZZ on the charge transfer ∆ Q . a bE σ ∗ AR σ ∗ ZZ electron accumulation while a negative value (blue) de-notes electron depletion. These bonding charge distri-butions clearly show the electron accumulation in themiddle region of M and X. The amount of charge local-ized in this region qualitatively indicates the strength ofthe covalent M-X bond. It is evident that as the com-position changes from sulphides to tellurides, the chargeaccumulation in the bonding regions gradually becomeless, indicating a weakening of M-X bond, lower Young’smodulus and ultimate strength.Fig. 4 also shows a large amount of charge transfer fromthe transition metal to the chalcogens. This charge redis-tribution can be quantitatively estimated by computingthe charge transfer from M to X. Table II presents themagnitude of the charge transfer obtained using Badercharge analysis. We depict the Young’s modulus E and ultimate strength σ ∗ as a function of charge trans-fer ∆ Q in Figs. 5(a) and 5(b), respectively. It is evi-dent that the mechanical properties of transition metaldichalcogenides exhibit a linearly increasing relation withthe charge transfer from transition metal atom to chalco-gen atoms. By fitting the values of the Young’s modulusand the ultimate strength with the charge transfer, thisrelationship can be simply described as E ( σ ∗ ) ∼ a + b ∆ Q where the values of fitted parameters a and b are listedin Table III. These expressions provide an approximateand simple description of the mechanical properties oftransition metal dichalcogenide monolayers sharing sim-ilar crystal structures. III. CONCLUSION
In summary, we have investigated the relationship be-tween the intrinsic mechanical response of the monolayer group VI TMDs to large elastic deformation, the atomic-level structure, and chemical composition. Our calcu-lations demonstrate that the chemical composition of amonolayer TMD strongly influences its mechanical re-sponse to large elastic deformation, with the chalcogensof W exhibiting much greater ultimate strength thanthe chalcogens of Mo. The stress response of monolayerTMDs also depends on the crystal symmetry, with thearmchair directions being consistently stronger across allchemical compositions. The origin of the observed me-chanical behavior can be attributed to a strong hybridiza-tion between the outermost p orbitals of the chalcogensand the d orbitals of the transition metal. This hybridiza-tion leads to a redistribution of the electronic chargeto the shared region between the transition metal andchalcogen atoms. Our study clearly highlights the inter-play between the atomic structure and the compositionof monolayer TMDs that lends them a great strengthto sustain large reversible deformations — deformationslarge enough to modify their electronic structure as en-visioned for novel optoelectronic and photonics applica-tions. It should be emphasized that our first-principlesinvestigation provides a realistic estimate of the intrinsicmechanical response of monolayer TMDs as measured bythe nano-indentation experiments. In such experiments,the stressed region of the monolayer under the atomicforce microscope tip can be expected to be defect-freeand, therefore, the measured stresses can approach theultimate strength. Therefore, our investigation of thestress-strain relation of monolayer TMDs presented herenot only yields crucial insights into mechanical behav-ior of these materials but also can be compared with theexperimental measurements. ACKNOWLEDGMENTS
J.L. and V.B.S. gratefully acknowledge the support ofthe Army Research Office through Contract W911NF-11-1-0171. N.V.M. acknowledges the computational supportfrom MASSIVE and NCI national facilities. ∗ [email protected] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang,Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A.Firsov, Science , 666 (2004). K. S. Novoselov, D. Jiang, F. Schedin, T. J. Booth, V. V.Khotkevich, S. V. Morozov, and A. K. Geim, Proc. Natl.Acad. Sci. , 10451 (2005). Y. Lin and J. W. Connell, Nanoscale , 6908 (2012). Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Coleman,and M. S. Strano, Nat. Nanotechnol. , 699 (2012). S. Balendhran, J. Deng, J. Z. Ou, S. Walia, J. Scott,J. Tang, K. L. Wang, M. R. Field, S. Russo, S. Zhuiykov,M. S. Strano, N. Medhekar, S. Sriram, M. Bhaskaran, andK. Kalantar-zadeh, Adv. Mater. , 109 (2013). J. N. Coleman, M. Lotya, A. O’Neill, S. D. Bergin, P. J.King, U. Khan, K. Young, A. Gaucher, S. De, R. J. Smith,I. V. Shvets, S. K. Arora, G. Stanton, H.-Y. Kim, K. Lee,G. T. Kim, G. S. Duesberg, T. Hallam, J. J. Boland,J. J. Wang, J. F. Donegan, J. C. Grunlan, G. Moriarty,A. Shmeliov, R. J. Nicholls, J. M. Perkins, E. M. Grieve- son, K. Theuwissen, D. W. McComb, P. D. Nellist, andV. Nicolosi, Science , 568 (2011). P. Joensen, R. Frindt, and S. Morrison, Mater. Res. Bull. , 457 (1986). A. Schumacher, L. Scandella, N. Kruse, and R. Prins,Surf. Sci. , L595 (1993). Y. Shi, W. Zhou, A.-Y. Lu, W. Fang, Y.-H. Lee, A. L. Hsu,S. M. Kim, K. K. Kim, H. Y. Yang, L.-J. Li, J.-C. Idrobo,and J. Kong, Nano Lett. , 2784 (2012). A. Castellanos-Gomez, M. Barkelid, S. Goossens, V. E.Calado, H. S. J. van der Zant, and G. A. Steele, NanoLett. , 3187 (2012). C. Mattevi, H. Kim, and M. Chhowalla, J. Mater. Chem. , 3324 (2011). J. Wilson and A. Yoffe, Adv. Phys. , 193 (1969). B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti,and A. Kis, Nat. Nanotechnol. , 147 (2011). A. Splendiani, L. Sun, Y. Zhang, T. Li, J. Kim, C.-Y.Chim, G. Galli, and F. Wang, Nano Lett. , 1271 (2010). G. Eda, H. Yamaguchi, D. Voiry, T. Fujita, M. Chen, andM. Chhowalla, Nano Lett. , 5111 (2011). J. K. Ellis, M. J. Lucero, and G. E. Scuseria, Appl. Phys.Lett. , 261908 (2011). Z. Y. Zhu, Y. C. Cheng, and U. Schwingenschl¨ogl, Phys.Rev. B , 153402 (2011). T. Cao, G. Wang, W. Han, H. Ye, C. Zhu, J. Shi, Q. Niu,P. Tan, E. Wang, B. Liu, and J. Feng, Nat. Comm. , 887(2012). D. Xiao, G.-B. Liu, W. Feng, X. Xu, and W. Yao, Phys.Rev. Lett. , 196802 (2012). H. Zeng, J. Dai, W. Yao, D. Xiao, and X. Cui, Nat. Nan-otechnol. (2012). K. F. Mak, K. He, J. Shan, and T. F. Heinz, Nat. Nan-otechnol. (2012). S. Bertolazzi, J. Brivio, and A. Kis, ACS Nano , 9703(2011). Q. Yue, J. Kang, Z. Shao, X. Zhang, S. Chang, G. Wang,S. Qin, and J. Li, Phys. Lett. A , 1166 (2012). W. S. Yun, S. W. Han, S. C. Hong, I. G. Kim, and J. D.Lee, Phys. Rev. B , 033305 (2012). P. Johari and V. B. Shenoy, ACS Nano , 5449 (2012). A. Castellanos-Gomez, M. Poot, G. A. Steele, H. S. J.van der Zant, N. Agra¨ıt, and G. Rubio-Bollinger, Adv.Mater. , 772 (2012). C. Ataca, M. Topsakal, E. Akturk, and S. Ciraci, J. Phys.Chem. C , 16354 (2011). X. Gonze, B. Amadon, P.-M. Anglade, J.-M. Beuken,F. Bottin, P. Boulanger, F. Bruneval, D. Caliste, R. Cara-cas, M. Cˆot´e, T. Deutsch, L. Genovese, P. Ghosez, M. Gi-antomassi, S. Goedecker, D. Hamann, P. Hermet, F. Jol-let, G. Jomard, S. Leroux, M. Mancini, S. Mazevet,M. Oliveira, G. Onida, Y. Pouillon, T. Rangel, G.-M. Rig-nanese, D. Sangalli, R. Shaltaf, M. Torrent, M. Verstraete,G. Zerah, and J. Zwanziger, Comput. Phys. Commun. , 2582 (2009). X. Gonze, G.-M. Rignanese, M. Verstraete, J.-M. Beuken,Y. Pouillon, R. Caracas, F. Jollet, M. Torrent, G. Zerah,M. Mikami, P. Ghosez, M. Veithen, J.-Y. Raty, V. Ole-vano, F. Bruneval, L. Reining, R. Godby, G. Onida, D. R.Hamann, and D. C. Allan, Zeit. Kristallogr. , 558(2005). C. Hartwigsen, S. Goedecker, and J. Hutter, Phys. Rev.B , 3641 (1998). A. Kumar and P. K. Ahluwalia, Eur. Phys. J. B , 186(2012). Y. Ding, Y. Wang, J. Ni, L. Shi, S. Shi, and W. Tang,Physica B: Condens. Matter , 2254 (2011). T. Lorenz, D. Teich, J.-O. Joswig, and G. Seifert, J. Phys.Chem. C , 11714 (2012). F. Liu, P. Ming, and J. Li, Phys. Rev. B , 064120 (2007). G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169(1996). P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994). R. F. W. Bader,
Atoms in Molecules: A Quantum Theory (Oxford University Press, USA, 1994). W. Tang, E. Sanville, and G. Henkelman, J. Phys.: Con-dens. Matter21