⇦Chinese Physics Letters, 2018, Vol. 35, No. 1, Article code 017302 Strain Effects on Properties of Phosphorene and Phosphorene Nanoribbons: a DFT and Tight Binding Study * Ruo-Yu Zhang(张若愚)1, Ji-Ming Zheng(郑继明)2**, Zhen-Yi Jiang(姜振益)1** Affiliations 1Department of Physics, Northwest University, Xi'an 710069 2National Key Laboratory of Photoelectric Technology and Functional Materials (Culture Base) in Shaanxi Province, National Photoelectric Technology and Functional Materials & Application of Science and Technology International Cooperation Base, Institute of Photonics & Photon-Technology, Northwest University, Xi'an 710069 Received 28 August 2017 *Supported by the National Natural Science Foundation of China under Grant Nos 51572219 and 11447030, the Natural Science Foundation of Shaanxi Province of China under Grant Nos 2014JM2-1008 and 2015JM1018, and the State Key Laboratory of Transient Optics and Photonics Technology 2015 Annual Open Fund under Grant No SKLST200915.
**Corresponding author. Email: zjm@nwu.edu.cn Citation Text: Zhang R Y, Zheng J M and Jiang Z Y 2018 Chin. Phys. Lett. 35 017302 Abstract We perform comprehensive density functional theory calculations of strain effect on electronic structure of black phosphorus (BP) and on BP nanoribbons. Both uniaxial and biaxial strain are applied, and the dramatic change of BP's band structure is observed. Under 0–8% uniaxial strain, the band gap can be modulated in the range of 0.55–1.06 eV, and a direct–indirect band gap transition causes strain over 4% in the $y$ direction. Under 0–8% biaxial strain, the band gap can be modulated in the range of 0.35–1.09 eV, and the band gap maintains directly. Applying strain to BP nanoribbon, the band gap value reduces or enlarges markedly either zigzag nanoribbon or armchair nanoribbon. Analyzing the orbital composition and using a tight-binding model we ascribe this band gap behavior to the competition between effects of different bond lengths on band gap. These results would enhance our understanding on strain effects on properties of BP and phosphorene nanoribbon. DOI:10.1088/0256-307X/35/1/017302 PACS:73.22.-f, 73.90.+f, 73.20.-r © 2018 Chinese Physics Society Article Text Two-dimensional (2D) materials, such as graphene and molybdenum disulfide (MoS$_{2}$), have shown potential application in biochemical sensors, electronic, and optoelectronic devices owing to their unique properties.[1-5] Very recently, a new 2D crystal has been successfully produced by mechanical exfoliation from bulk black phosphorus (BP), and named as phosphorene.[6,7] Similar to transition metal dichalcogenides (TMDs), BP also exhibits a layer-dependent band gap whose range is 0.3 eV in bulk and 2 eV in monolayer.[8,9] However, different from TMDs, phosphorene has much larger carrier mobility (around 1000 cm$^{2}$/V$\cdot$s, while it is around only 200 cm$^{2}$/V$\cdot$s in MoS$_{2}$), and shows large anisotropic mobility of charge carriers, and optic properties.[10,11] The features indicate that phosphorene is a promising candidate for the next generation of electronic and optical applications. To make phosphorene more multifunctional, it is important to explore efficient ways to modulate its properties. There have been plenty of works on the band structure tuning by mechanical stress,[12] surface decoration,[13] electric field,[14] doping,[15] and so on. Stress has been used for a long time to tune electronic structure of semiconductors, and has also been proved to be an effective and simple means of modulating band gap in 2D materials.[16-19] In particular, 2D crystals can sustain much larger strains than their bulk crystals, for example, monolayer graphene and MoS$_{2}$ have been stretched to their intrinsic limits ($\sim$15% for graphene and $\sim$11% for MoS$_2$)[20-22] without substantial damage. It has also been reported that phosphorene can sustain a large strain to 30% without lattice deformation and that its band structures are also extremely sensitive to strain.[19] Stress is inevitable. In addition to a mechanical way, cutting BP into ribbons, stacking BP on other 2D materials to form van der Waals (vdW) heterostructure, applying external fields can also induce the stress effect on BP when integrating black phosphorus in the device. Nevertheless, this is not just a defect. Researchers have managed to understand it and use this effect. Peng et al.[19] have revealed the origin of the gap transition under strain according to the bond nature of monolayer BP. Using density functional theory (DFT) and tight-binding models, Rodin et al.[23] predicted a semiconductor-to-metal transition in phosphorene under compression. Yang's group demonstrated that unique anisotropic free-carrier mobility of BP can be controlled by using simple strain conditions. With the appropriate biaxial or uniaxial strain (4–6%), they rotate the preferred conducting direction by 90$^{\circ}$. With what is commonly used in graphene-based devices, cutting BP into ribbons will be a routine to reduce the BP-based device scale. Thus quantum confinement and edge effects will play crucial roles in determining the band gap,[24] and hence may induce different band gap behaviors from that observed in bulk states under stress. However, there is still a lack of studies about the strain effects of the BP ribbons.[25,26] In this work, first principle methods are used to systematically study the strain effects on the phosphorene electronic structure, and on the corresponding nanoribbons. We find that the band gap of phosphorene increases (or decreases) when the stress is increased (or decreased) with uniaxial or biaxial strain, and a direct to indirect band gap transition is also observed. Lastly we discuss the reason for the change of band gap and give a brief conclusion. The optimized geometrical structure, and electronic properties of BP and its nanoribbons all are exploited using density functional theory (DFT) as implemented in the Vienna ab initio package simulation code (VASP).[27,28] The electron exchange and correlation are described by the generalized gradient approximation of Perdew–Burke–Ernzerhof (PBE) functions.[29] The convergence tests of the $k$-points sampling and energy cutoff have been carefully examined before the calculation of the electronic structure. The plane-wave energy cutoff is set to be 400 eV and the $k$-point sampling grid is $14\times10\times1$ in monolayer phosphorus, and $31\times1\times1$ in monolayer phosphorus nanoribbons, respectively. Vacuum spacing between neighboring layers is set to be 20 Å to avoid the interaction between nanoribbons. The total ground state energy is converged within 10$^{-6}$ eV per unit cell. Monolayer phosphorus is relaxed and the final force exerted on each P atom is less than 0.01 eV/Å. Here 30 points are taken along each high symmetry line for calculating band structure.
cpl-35-1-017302-fig1.png
Fig. 1. The ball-stick models [(a), (b)] of monolayer BP. The dashed rectangle in (b) indicates a unit cell. (c) DFT-GGA calculated band structure of BP. The Fermi level is set to 0 eV.
Structure of BP is shown in Figs. 1(a) and 1(b). The pristine BP is built with a $5\times5$ supercell consistency, which contains 100 atoms in total. Unlike a flat structure of many 2D materials, single layer BP has a puckered honeycomb structure with each P atom covalently bonded with adjacent P atoms. Lattice constants for phosphorene after relaxing are $a=3.30$ Å and $b=4.61$ Å. We define the applied strain as $\varepsilon _x =\frac{a_x -a}{a}$, $\varepsilon _y =\frac{b_y -b}{b}$, where $a_{x}$ and $b_{y}$ are the lattice constants along the $x$ ($y$) direction after applying tensile. The lattice constant along the direction with strain is fully relaxed through both energy and force minimization to ensure a stable structure. In our DFT calculations, monolayer BP is an approximately direct-gap two-dimensional material (the energy of actually highest point within the conduction band is merely 0.1 eV higher than point ${\it \Gamma}$). As shown in Fig. 1(c), there is a 0.92 eV band gap in BP, which is in good agreement with the previous works. One of the most fantastic characteristics that can be observed in the BP band spectrum is the highly anisotropic band along $X$–${\it \Gamma}$–$Y$ direction in reciprocal space. In particular, the band dispersion along $Y$–${\it \Gamma}$ direction (armchair direction) is much steeper than that along $X$–${\it \Gamma}$ direction (zigzag direction), indicating the anisotropy of BP, such as effective mass which is much higher in the armchair direction (0.15$m_{\rm e}$, where $m_{\rm e}$ is the bare electron mass) than that in the zigzag direction (1.25$m_{\rm e}$). Figure 2 shows the band structure of phosphorene under an increasing uniaxial strain along $x$ or $y$ direction, respectively. The most remarkable variation in these band structures occurs along $X$–${\it \Gamma}$ and ${\it \Gamma}$–$Y$ directions, especially in three points, $X (0.5, 0, 0)$, $Y (0, 0.5, 0)$ and ${\it \Gamma}$. Under strain increasing from 0% to 8% along the $X$ direction, there are two steps of band gap change, 0%–3% and 3%–8%. In the former step, the band gap is proportional to the strain (from initial 0.92 eV to the maximum value 1.06 eV), while in the latter, it is dropped very quickly (from 1.06 eV to 0.55 eV) as shown in Fig. 2(a). This phenomenon is originated from the downward shift of the conduction band minimum (CBM), and upward shift of valence band maximum (VBM). It is worth noting that monolayer BP is always a direct band gap semiconductor in these processes.
cpl-35-1-017302-fig2.png
Fig. 2. Calculated band structures of monolayer BP (a) for applied strain along $x$ direction and (b) for applied strain along $y$ direction from 0% to 8%.
However, when the stress is applied along the $y$ direction, the result is completely different. There is a direct to indirect band gap transition in Fig. 2(b). CBM is at ${\it \Gamma}$ point when strain is less than 4% (the energies of A and B are 0.51 eV and 0.55 eV, respectively), while it transfers from ${\it \Gamma}$ to $Y$ when the strain is 4 % (the energies of C and D are 0.53 eV and 0.51 eV). We could believe that the energies of two points competing for CMB are equal for strain 3%–4%. Thus point D becomes CBM for strain $>$4%. In the direct band gap stage, the value of the band gap grows from 0.92 eV to 1.03 eV while in the indirect band gap stage it only grows to 0.93 eV. In this process, the valence band almost does not change and the band gap value depends on CBM. When biaxial strain is added to the monolayer BP, Fig. 3(a) shows the BP's band structures when it is under biaxial strain range from $-$5% to +5%. The data of band gap change are also given in Fig. 3(b). The lowest band gap value 0.34 eV appears at the ${\it \Gamma}$ point for $-$5% strain, while it increases to the maximum value 1.09 eV for +2% stress applied. With the increasing strain to +5%, the band gap reaches 0.84 eV. In the whole process of increasing biaxial stress, the monolayer BP is a direct band gap semiconductor and both CBM and VBM are located at ${\it \Gamma}$. VBM is reduced from $-$5%–2% then increases from 2%–5%, while CBM is opposite.
cpl-35-1-017302-fig3.png
Fig. 3. (a) Calculated band structures of phosphorene with increasing biaxial strain. (b) Band gap of phosphorene as a function of the uniaxial strain $\varepsilon_{x}$ and $\varepsilon_{y}$, respectively.
cpl-35-1-017302-fig4.png
Fig. 4. Models of pristine BP nanoribbons along (a) armchair and (b) zigzag orientations. Side views of pristine BP nanoribbons along (c) armchair and (d) zigzag orientations. Both of them are hydrogenated at the edge.
By cutting phosphorene along two vertical orientations, respectively, armchair or zigzag nanoribbons are formed as shown in Fig. 4. The strain effects on electronic structure of these ribbons are summarized in Fig. 5. In contrast to pristine phosphorene which is a direct band gap semiconductor, phosphorene nanoribbon with bare edges exhibits metallic properties due to the dangling edge states. These edge states are not stable and should be saturated with hydrogen atoms. After geometrical relaxation, hydrogen saturation leads the ribbon to open a direct band gap, 1.298 eV and 1.035 eV for armchair or zigzag nanoribbons, respectively. Both of them are direct band gap semiconductors, with VBM and CBM located at the ${\it \Gamma}$ point. As shown in Fig. 5, BP nanoribbons maintain direct band gap at the ${\it \Gamma}$ point under all strained states. For zigzag nanoribbons, in Fig. 5(a), the values of the band gap decrease monotonically to the minimum value of 0.28 eV at $x=10{\%}$, as the strain increases. However, for armchair nanoribbons, the band gap values change differently as the strain increases. Two steps can be observed in Fig. 5(b). At the first step (when strain $x\leq8{\%}$), the band gaps increase from 1.04 eV to 1.22 eV, at the second step (when strain $8{\%}\leq x\leq10{\%}$), they decrease down to 1.18 eV at $x=10{\%}$. Energy of the highest energy valance band (redline) is almost unchanged throughout the process of stress. As for the conduction band, there are two points competing for CBM. The first conduction band is steep around the ${\it \Gamma}$ point (blue line when $x\leq8{\%}$), while the second one (blue line when $8{\%}\leq x\leq10{\%}$) is much more flattened. CBM is located from the steeper to the flatter conduction band. Energy of the first conduction band rises with the increasing positive strain while the energy of second conduction band hardly changes. At 8% strain, CBM is located from the first conduction band to the second one.
cpl-35-1-017302-fig5.png
Fig. 5. Calculated band structures of monolayer BP nanoribbons for (a) zigzag nanoribbon and (b) armchair nanoribbon. Here we apply the strain from 0% to 10%.
Energy bands of materials are composed of different partitions of atomic orbital and these orbitals respond differently to strains along different directions. To understand the band shift when phosphorene nanoribbons are under stress, it is helpful to analyze the behavior of these orbitals under strains. Figure 6 shows the PDOS of phosphorene nanoribbons. The VBM is mainly composed of $p_{z}$ orbitals, while the CBM shows mixed features between $p_{z}$ and $p_{y}$ orbitals. As the stress on phosphorene increases, its band structure gradually changes to be indirect, and CBM now shows a mixture of $p_{x}$ and $s$ orbitals. In zigzag nanoribbons, both CBM and VBM have strong contributions from $p_{z}$ orbital, while in armchair nanoribbons, we concentrate on VBM and two lowest energy points in the first and second conduction bands. The VBM and CBM when the nanoribbons are applied tensile less than $x=8{\%}$ consist of $p_{z}$ orbitals and the lowest energy point of second conduction band is made up of $p_{y}$ orbital.
cpl-35-1-017302-fig6.png
Fig. 6. Partial density of states (PDOS) in BP nanoribbons for (a) zigzag nanoribbon and (b) armchair nanoribbon.
To explore the critical strain in which the increasing–decreasing band gap value transition occurs, we exploit general mechanism to explain the CBM and VBM converting and varying by means of Heitler–London's exchange energy model. Shifting of energy bands has been well illustrated in monolayer phosphorus. In this model, energy changing of the bonding and anti-bonding states could be described by $$ E=2E_0+\Delta E, $$ where $E_{0}$ is the ground-state energy of single phosphorus atom, and $\Delta E$ is the energy change due to the relative position change between the two phosphorus atoms. There are two possibilities of $\Delta E$ depending upon the bonding types. The energy is $E=\frac{C(R)\mp X(R)}{1-S(R)}$ for the bonding type and anti-bonding type, respectively, where $S$, $C$ and $X$ all are functions of atomic distance, called the overlap integral, the Coulomb integral and the exchange integral, respectively. These integrals have a simple interpretation: $S$ is just the overlap between neighbor phosphorus atoms of the two wave functions which is much smaller than 1, thus we could ignore it in the equation. Moreover, according to the formula, bonding energy also increases and decreases for anti-bonding energy as $U$ increases, while non-bonding energy hardly changes. Now we will take advantage of this rule to explain band gap changes in nanoribbon. As shown in Fig. 7, in armchair nanoribbon, orbitals of VBM are non-bonding while in CBM they are anti-bonding. In zigzag nanoribbon, orbitals of VBM are bonding while in CBM they are anti-bonding. The bond nature of these states could fit well with bonding energy law according to the formula given above.
cpl-35-1-017302-fig7.png
Fig. 7. The electron density contour plots of (a) VBM of zigzag phosphorene nanoribbon, (b) CBM of zigzag phosphorene nanoribbon, (c) VBM of armchair phosphorene nanoribbon, and (d) CBM (0%–8%) of armchair phosphorene nanoribbon.
The tight-binding Hamiltonian for phosphorene has been proposed by Rudenko and Katsnelson.[30] The model is $H=\sum\nolimits_{i,j} {t_{i,j}} c_i^† c_j$, where $c_i^†$ and $c_j$ represent the creation of one electron in the $i$th lattice point and the annihilation of one electron in the $j$th lattice point, $t_{i,j}$ are the hopping integrals between $i$th and $j$th lattice points, the summation is over all the lattice points of the ribbon, and the hopping integrals in the above Hamiltonian are explicitly taken as $t_{1}=-1.220$ eV, $t_{2}=3.665$ eV, $t_{3}=-0.205$ eV, $t_{4}=-0.105$ eV, and $t_{5}=-0.055$ eV, respectively. Tight-binding model could give a better understanding by capturing the information around the ${\it \Gamma}$ point.[31] When the ribbons are under strain states, atoms are gradually relaxed to new site and bond length changes correspondingly, which also lead to the variation of the hopping integrals $t_{i,j}$.[32] Among all the hopping integrals, $t_{1}$ and $t_{2}$ are the most important ones, which are the hopping integral along zigzag chains ($t_{1}$) and between the nearest atoms of the different layers ($t_{2}$). The change of the two parameters has different effects on the band gap of the ribbon. Figure 8(a) shows the details in armchair and zigzag BPRs, where 'change of hopping integrals' means deceasing percentage of $t_{1}$ or $t_{2}$ compared with the original ones $(\frac{\Delta t_{1/2} }{t_{1/2} }\times 100\% )$. It is clear that the band gap increases (decreases) with decreasing $t_{1}$ ($t_{2}$). With tensile strain increasing, the bond length increases correspondingly. Figure 8(b) gives the variation of $t_{1}/t_{2}$ under the stretching. The bond length corresponding to $t_{1}$ increases quickly at first stretching, then keeps nearly unchanged, and becomes to shrink when the strain is larger than 8%. The bond length corresponding to $t_{2}$ increases linearly under all the strain states. Comparing this with our DFT results, the band gap variation would be ascribed to competition between influence of $t_{1}$ and $t_{2}$ on the band gap. At first step of stretching, $t_{1}$ is the dominant parameter which leads to the increase of the band gap with increasing the strain. When the strain is larger than 8%, $t_{2}$ becomes the controlling factor and leads to the decrease of the band gap.
cpl-35-1-017302-fig8.png
Fig. 8. (a) Band gap under variation of $t_{1}$ (left) and $t_{2}$ (right) as a function of change of hopping integrals. (b) Increase of bond length corresponding to $t_{1}$ (left) and $t_{2}$ (right) as a function of change of hopping integrals.
In summary, we have investigated the influence of uniaxial and biaxial strains on electronic structures of monolayer BP and monolayer passivate BP nanoribbon. The band gaps of monolayer BP and BP nanoribbon dramatically change with the axial strain. The band gap of monolayer phosphorene has direct–indirect–direct transitions with uniaxial strain. In BP nanoribbon, it has a competition between two conduction bands during applying stress. We develop a general mechanism to explain the different energy shifts with strain based on Heitler–London's exchange energy model. It is found that energy change, which is related to strain, is closely related to the bonding/reaction properties of the orbit. The increase or decrease of atomic P–P distance (corresponding to a positive or negative tensile strain) causes the bonding energy $E_{\rm bonding}$ to increase and the anti-bonding energy $E_{\rm anti-bonding}$ to decrease. The bond nature of these states can be used to explain their energy variation. We have proved the mechanism of bonding energy to explain the process according to the bond nature of their electronic orbitals in BP nanoribbon and tight bonding models. BP has shown great potential in optical and electronic fields due to its direct band gap. We demonstrate that the band gap value can be tuned by strain in BP or BP nanoribbon and even achieve a indirect–direct band gap process. References Electric Field Effect in Atomically Thin Carbon FilmsControl of Graphene's Properties by Reversible Hydrogenation: Evidence for GraphaneEmerging Photoluminescence in Monolayer MoS 2Vapour phase growth and grain boundary structure of molybdenum disulphide atomic layersElectronic Muscles and Skins: A Review of Soft Sensors and ActuatorsPhosphorene: An Unexplored 2D Semiconductor with a High Hole MobilityBlack phosphorus field-effect transistorsAb initio studies on atomic and electronic structures of black phosphorusHigh-mobility transport anisotropy and linear dichroism in few-layer black phosphorusPhosphorene Nanoribbons, Phosphorus Nanotubes, and van der Waals MultilayersMeasurement of mobility in dual-gated MoS2 transistorsStrain and Orientation Modulated Bandgaps and Effective Masses of Phosphorene NanoribbonsA First-Principles Study on Electron Donor and Acceptor Molecules Adsorbed on PhosphoreneStrain-Engineering the Anisotropic Electrical Conductance of Few-Layer Black PhosphorusTransition Metal Doped Phosphorene: First-Principles StudyTunable electronic structure in stained two dimensional van der Waals g-C 2 N/ X Se 2 ( X = Mo, W) heterostructuresStrain-induced semiconductor to metal transition in few-layer black phosphorus from first principlesBand engineering of dichalcogenide MX2 nanosheets (M = Mo, W and X = S, Se) by out-of-plane pressureStrain-engineered direct-indirect band gap transition and its mechanism in two-dimensional phosphoreneLarge-scale pattern growth of graphene films for stretchable transparent electrodesMeasurement of the Elastic Properties and Intrinsic Strength of Monolayer GrapheneMechanical properties of freely suspended semiconducting graphene-like layers based on MoS2Strain-Induced Gap Modification in Black PhosphorusEnergy Gaps in Graphene NanoribbonsAb initio studies of thermodynamic and electronic properties of phosphorene nanoribbonsElectronic Properties of Edge-Hydrogenated Phosphorene Nanoribbons: A First-Principles StudyEfficient iterative schemes for ab initio total-energy calculations using a plane-wave basis setEfficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis setGeneralized Gradient Approximation Made SimpleQuasiparticle band structure and tight-binding model for single- and bilayer black phosphorusMultiple unpinned Dirac points in group-Va single-layers with phosphorene structureStrain modification on electronic transport of the phosphorene nanoribbon
[1] Novoselov K S, Geim A K, Morozov S V, Jiang D, Zhang Y, Dubonos S V, Grigorieva I V and Firsov A A 2004 Science 306 666
[2] Elias D C, Nair R R, Mohiuddin T M, Morozov S V, Blake P, Halsall M P, Ferrari A C, Boukhvalov D W, Katsnelson M I, Geim A K and Novoselov K S 2009 Science 323 610
[3] Splendiani A, Sun L, Zhang Y, Li T, Kim J, Chim C Y, Galli G and Wang F 2010 Nano Lett. 10 1271
[4] Najmaei S, Liu Z, Zhou W, Zou X, Shi G, Lei S, Yakobson B I, Idrobo J C, Ajayan P M and Lou J 2013 Nat. Mater. 12 754
[5] Chen D and Pei Q 2017 Chem. Rev. 117 11239
[6] Liu H, Neal A T, Zhu Z, Luo Z, Xu X, Tománek D and Ye P D 2014 ACS Nano 8 4033
[7] Li L, Yu Y, Ye G J, Ge Q, Ou X, Wu H, Feng D, Chen X H and Zhang Y 2014 Nat. Nanotechnol. 9 372
[8] Du Y, Ouyang C, Shi S and Lei M 2010 J. Appl. Phys. 107 093718
[9] Qiao J, Kong X, Hu Z X, Yang F and Ji W 2014 Nat. Commun. 5 4475
[10] Guo H, Lu N, Dai J, Wu X and Zeng X C 2014 J. Phys. Chem. C 118 14051
[11] Fuhrer M S and Hone J 2013 Nat. Nanotechnol. 8 146
[12] Han X, Morgan Stewart H, Shevlin S A, Catlow C R A and Guo Z X 2014 Nano Lett. 14 4607
[13] Zhang R, Li B and Yang J 2015 J. Phys. Chem. C 119 2871
[14] Fei R Y L 2014 Nano Lett. 14 2884
[15] Hashmi A and Hong J 2015 J. Phys. Chem. C 119 9198
[16] Zheng Z D, Wang X C and Mi W B 2017 Physica E 94 148
[17] Ju W, Li T, Wang H, Yong Y and Sun J 2015 Chem. Phys. Lett. 622 109
[18] Su X, Zhang R, Guo C, Zheng J and Ren Z 2014 Phys. Lett. A 378 745
[19] Peng X, Wei Q and Copple A 2014 Phys. Rev. B 90 085402
[20] Kim K S, Zhao Y, Jang H, Lee S Y, Kim J M, Kim K S, Ahn J H, Kim P, Choi J Y and Hong B H 2009 Nature 457 706
[21] Lee C, Wei X, Kysar J W and Hone J 2008 Science 321 385
[22] Castellanos Gomez A, Poot M, Steele G A, Zant H S J V D, Agraït N and Rubio Bollinger G 2012 Nanoscale Res. Lett. 7 233
[23] Rodin A S, Carvalho A and Castro Neto A H 2014 Phys. Rev. Lett. 112 176801
[24] Son Y W, Cohen M L and Louie S G 2006 Phys. Rev. Lett. 97 216803
[25] Ramasubramaniam A and Muniz A R 2014 Phys. Rev. B 90 085424
[26] Li W, Zhang G and Zhang Y W 2014 J. Phys. Chem. C 118 22368
[27] Kresse G and Furthmüller J 1996 Phys. Rev. B 54 11169
[28] Kresse G and Furthmüller J 1996 Comput. Mater. Sci. 6 15
[29] Perdew J P, Burke K and Ernzerhof M 1996 Phys. Rev. Lett. 77 3865
[30] Rudenko A N and Katsnelson M I 2014 Phys. Rev. B 89 201408
[31] Lu Y, Zhou D, Chang G, Guan S, Chen W, Jiang Y, Jiang J, Wang X s, Yang S A, Feng Y P, Kawazoe Y and Lin H 2016 NPJ Comput. Mater. 2 16011
[32] Yuan Y W and Cheng F 2017 AIP Adv. 7 075310