Robust quasi-ohmic contact against angle rotation in noble transition-metal-dichalcogenide/graphene heterobilayers

Abstract

Stacked 2D crystals provide great possibilities to achieve desired functionalities with varying thickness and composition properties. In this work, the electronic properties of heterobilayers of noble transition-metal-dichalcogenide (PdS₂, PtS₂)/graphene with interlayer rotation angles are systematically investigated by density functional theory calculations. The two layers are combined by weak van der Waals interactions. Each component preserves its band structure in the hybrid system. A charge transference occurs at the interface of the bilayers leading to p-type doping for graphene and n-type doping for PdS₂ (PtS₂). Schottky barriers ranging from 0 eV to 0.057 eV appear in both heterobilayers with different orientations between the two monolayers. These results suggest that PdS₂ (PtS₂)/graphene heterobilayers have quasi-ohmic contact against interface angle rotation.

---

1. Introduction

Since the successful invention of isolated graphene in experiments, two-dimensional (2D) materials have received much attention.^1–5 Layered transition-metal dichalcogenide (TMDC) nanostructures, which consist of one layer of transition metal atoms (Mo, W, V, or Ti) sandwiched between two layers of chalcogen atoms (S, Se, or Te), are now rapidly emerging, as they display unique chemical and physical properties that are absent or difficult to obtain in other 2D materials.⁶ As many TMDCs are inherent semiconductors by nature, a wide range of TMDC based applications in electronics, optoelectronics, photovoltaics, and catalysis have been demonstrated.^7–13 While intensive interests are focused on MoS₂ and WSe₂, a new type of TMDC, noble-transition-metal dichalcogenides (nTMDCs), remains rarely explored, and understanding the physics behind this system is fundamental to their practical applications. The pioneering work by Grønvold et al. indicates that the noble metals Pt and Pd can form layered structures with S atoms.^14–16 These monolayer noble-transition-metal disulfides MS₂ (M = Pt and Pd) are typical semiconductors with indirect band gaps of 1.75 eV for PtS₂ and 1.11 eV for PdS₂.¹⁷ The PtS₂ crystal has a layer-dependent indirect bandgap from 1.6 eV (monolayer) to 0.25 eV (bulk).¹⁸ Moreover, the phonon-limited mobility of platinum disulfide (PtS₂) can be as high as 1107 cm² V⁻¹ s⁻¹ at room temperature, about four times the carrier mobility of MoS₂.¹⁹ Therefore, PtS₂ may show better performance than MoS₂ when it serves as a channel in field effect transistors (FETs).

Recent efforts have focused on taking advantage of the individual properties of different 2D materials by fabricating heterostructures, which are vertical stacks of 2D layers of dissimilar materials held together by van der Waals (vdW) forces.²⁰ Graphene is an ideal component of multilayer heterojunctions including graphene/silicene, graphene/hexagonal boron nitride, graphene/stanene, graphene/phosphorene, and graphene/MoS₂,^21–26 due to its high stability and outstanding electronic properties. Such heterostructures may lead to the opening of a bandgap in graphene without impairing its electronic mobility.^27–29 Moreover, graphene is a suitable source electrode for TMDC based FETs. PtSe₂/G has a tunable layer-dependent Schottky barrier.³⁰ The traditional backgate MoS₂ field effect transistor with ion gel dielectrics has a rather low carrier mobility of 0.5–20 cm² V⁻¹ s⁻¹, while the characteristics of high mobility from graphene and high on-off ratio from MoS₂ are properly balanced in the MoS₂/graphene heterojunction based FET with an on–off ratio of 105 and a carrier mobility of 58.7 cm² V⁻¹ s⁻¹.^31–33 Different interlayer orientations might form in the synthesis of multilayer 2D materials, and interlayer orientation may affect the electronic properties.^34–41 To the best of our knowledge, the influence of interlayer rotation on the electronic properties of nTMDC/graphene hybrids is unknown.

In this paper, we have investigated the structure and electronic properties of nTMDC/graphene heterojunctions (PdS₂/G and PtS₂/G). Furthermore, we study the effect of interlayer orientation between two different 2D materials as a tuning parameter for the electronic properties of the bilayer heterostructure. The relative orientation of graphene with respect to the MS₂ layer apparently influences the interlayer spacing and the binding energy, but only slightly influences the value of the bandgap in MS₂. Robust quasi-ohmic contact against angle rotation is found in both PtS₂/G and PdS₂/G heterostructures.

2. Methods and computational details

Our calculations were performed using the Vienna ab initio simulation package (VASP) within the framework of the plane-wave pseudopotential DFT method.^42,43 The ion–electron interactions were described by the projected augmented wave method (PAW), and the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE) was employed to describe the exchange–correlation interactions.^44–47 In addition, the van der Waals interactions were considered using the DFT-D2 method.⁴⁸ The Brillouin zone was sampled by a k-point mesh of 0.05 Å⁻¹ separation in reciprocal space within the Monkhorst–Pack scheme, and the kinetic energy cutoff was chosen to be 400 eV.⁴⁹ A vacuum region larger than 15 Å was added to avoid spurious interactions between periodic images. All atomic positions were fully optimized until the forces acting on each atom were less than 0.01 eV Å⁻¹.

Twisted bilayers were modeled using accidental angular commensurations.⁵⁰ In a hexagonal lattice whose basis vector is a₁ and a₂, a skewed supercell with a basis vector (na₁ + ma₂) has a corresponding skewed angle, . The MS₂/G hybrid sheets were designed by looking for pairs of commensurate skewed cells of the two lattices with similar periodicity. The rotation angle of graphene relative to PdS₂ (PtS₂) can be expressed as θ_G − θ_MS2, namely the angle between the two zigzag directions, as shown in Fig. 1. Heterostructures are modeled by expanding the two monolayer unit cells to search the proper lattice combinations. To reduce the impact of external strain on MS₂ as much as possible, we chose the MS′₂ lattice parameter as the lattice parameter of the studied system when modeling the heterostructures. To avoid a large amount of computation and to sustain the original stability of the atom layers, we searched the combinations which could reach the balance between a little mismatch and a small number of atoms. We restricted the mismatch range from −5% to 5%. In this way, we found six twisted superstructures with relatively smaller numbers of atoms, namely, . In these superstructures, the number of atoms inside the unit cell increases as the orientation angle increases. In order to study the difference between the zigzag (ZZ) direction and the armchair direction (AC), we need a superstructure where the AC direction of graphene coincides with the ZZ direction of the MS₂ layer. However, to reduce the computational complexity, we adopted four typical superstructures as shown in Fig. 1, namely, , corresponding to rotation angles of 0, 7.3, 13.9 and 30, respectively. The interlayer displacement is confined to a stacking order where a C atom of graphene and a S atom of MS₂ are superimposed within each unit cell, considering that the interlayer bonding is insensitive to their mutual translation.

Fig. 1 Atomic geometries of the MS₂/G heterobilayer (M = Pt, Pd). (a and c–e) Top views of MS₂/G with rotation angles of 0, 7.3, 13.9 and 30 degrees, respectively. (b) Side view of the MS₂/G heterojunction, Δd indicates the interlayer distance between the two monolayers. The yellow, blue and silver spheres indicate sulfur, M and carbon atoms, respectively.

3. Results and discussions

3.1 Structure and energetics

The key geometrical properties of the two 1 × 1 non-twisted configurations are listed in Table 1. For PtS₂/G, the lowest S atom layer lies about 3.51 Å below the graphene. For PdS₂/G, the interface distance is about 3.45 Å. These interlayer distance values are close to those of the van der Waals systems such as graphite. Therefore, we infer that the interactions between graphene and PtS₂ or PdS₂ in the heterostructures are relatively weak. The interlayer distance has little change when the orientation is considered in a range of 0.01 Å, and this is very similar to twisted MoS₂/graphene heterobilayers.⁴¹ To quantitatively evaluate the interlayer interaction in the heterojunction, we calculated the binding energy (E_b), based on the following equation:

E_b = (E_G + E_MS2 − E_MS2/G)/N_C

where E_G is the total energy of graphene, E_MS2 is the total energy of MS₂, E_MS2/G is the energy of the MS₂/G heterobilayer, and N_C is the number of carbon atoms. The average binding energies are about 265 meV for PtS₂/G and 245 meV for PdS₂/G per C atom, respectively. Therefore, the PtS₂/G and PdS₂/G heterostructures are stable against thermal fluctuations at room temperature. The tendency of calculated E_b is plotted in Fig. 2. Interface interactions become stronger with the increase of the angles. When the rotation angle increases to 30°, the BE reaches a maximum value of 308 (248) meV for PtS₂/G (PdS₂/G), similar to the trend of the BE in the MoS₂/G system, which increases from 199 meV to 226 meV when the angle increases from 3 degrees to 25.9 degrees.⁴¹ Therefore, the system is more stable when the ZZ direction of graphene coincides with the AC direction of MS₂. The effect of the rotation angle on the graphene’s thickness and bond length is also almost negligible.

Table 1 Calculated properties of the two kinds of heterojunction. a denotes the lattice parameter of certain primitive cells; Δd denotes the surface distance between two layers; h denotes the height of the MS₂ layer; Φ_TB denotes the tunneling barrier of certain heterojunctions

	θ (°)	a (Å)	Mismatch	Δd (Å)	h (Å)	ΦTB (eV)
PtS2/G	0	7.19	−2.3%	3.51	2.65	3.95
	7.3	15.66	1.9%	3.51	2.72	4.23
	13.9	12.95	1.4%	3.55	2.66	3.75
	30	17.96	5.2%	3.51	2.66	4.40
PdS2/G	0	7.07	−4.5%	3.45	2.74	3.64
	7.3	15.40	1.9%	3.52	2.74	3.29
	13.9	12.74	1.4%	3.53	2.75	3.48
	30	17.66	3.6%	3.55	2.75	2.50

---

Fig. 2 Variation of binding energy of stable PdS₂/G (PtS₂/G) heterojunctions with the interface angle.

3.2 Electronic properties

The electronic properties of graphene-related heterostructures might be affected by interlayer rotations. To reveal this effect on MS₂, the band structures of the heterostructures with different rotation angles have been calculated as shown in Fig. 3. The PtS₂ and PdS₂ band gaps in the heterostructures will be greater than the current PBE values in HSE, while the band structure of graphene will not change significantly.⁵⁷ The spin orbital effect is not considered because no significant changes in the band gap or the band structure were observed in any of the monolayers of PtS₂ and PdS₂, and the band gap has no obvious change and the positions of the CBM and VBM do not move.¹⁷ To distinguish the contribution of each constituent layer more clearly, the energy states originated primarily from graphene are highlighted by red lines. All the Dirac points of graphene move to gamma points because the primitive cells of graphene used are multiples of . The free-standing PtS₂ and PdS₂ monolayers are semiconductors with bandgaps of 1.11 eV and 1.75 eV, respectively.¹⁷ When in contact with graphene, the band structure of PtS₂ (PdS₂) and graphene can be easily recognized, indicating the weak interlayer interaction. When the interlayer rotation angle is zero, the band gaps of PtS₂ and PdS₂ are 0.95 eV and 1.71 eV, respectively. The reduction of the bandgap of PtS₂ (PdS₂) after contact with the graphene layer is caused by the weak interlayer interactions. Remarkably, the Fermi level of the heterojunction moves to the position near the conduct band minimum (CBM) of MS₂, indicating that MS₂ is n-doped. The DFT calculations also show that MoS₂ and WS₂ have been slightly n-type doped in the MoS₂/G and WS₂/G heterostructures.^51,52 In the experimentally fabricated MoS₂ and WS₂ transistors with graphene electrodes, such n-type transfer characteristics have also been observed.^53,54 We believe this theoretical and experimental consistency with similar systems to our studied ones provides strong support for the reliability of our prediction. Now we can check the effect of different interlayer rotation angles. In both PtS₂/G and PdS₂/G heterostrucures, graphene always retains its semi-metallic character, and there is no band gap opened at its Dirac point, as the interlayer rotation angle changes. An obvious variation of the MS₂ band structure is the positions of the valence band maximum (VBM) and conduct band minimum (CBM). When the interlayer angle is greater than 7 degrees, the CBM moves towards the K point in PdS₂/G. For PtS₂/G, the CBM and VBM all locate between the K and gamma points, except for the interlayer rotation angles, which are 13.9 degrees. As the rotation angle increases, the band gap of PdS₂ is almost unchanged while that of PtS₂ changes significantly from 1.71 eV to 1.56 eV when in contact with graphene.

Fig. 3 Band structure of the monolayers PdS₂ (a–d) and PtS₂ (e–h) on graphene with different interface angles. The red lines and blue lines represent the band structures of graphene and MS₂, respectively. The Fermi level is set to zero and marked by dashed lines.

Due to the quasi-metal nature of graphene, this kind of heterojunction can be regarded as a metal–semiconductor contact. For n-type semiconductors, electrons near the CBM are primarily responsible for electric conduction. At the metal–semiconductor interface, the electrons of lowest-lying states (located at the CBM) for an n-type semiconductor can spontaneously flow to the metal. The states will be offset when electronic transport occurs across the metal–semiconductor interface known as the n-type Schottky barrier height (SBH). The flow of electrons from the semiconductor to the metal is easier than conduction in the opposite direction. The SBH controls the transport of electrons crossing the metal–semiconductor interface. Because the current flow across the MS₂ interface depends on the magnitude of the SBH, the SBH is clearly the most important property of a metal–semiconductor interface. The vertical Schottky barrier Φ_V for the above case can be obtained from the energy difference between the E_f value of the interfacial system and the CBM (electron SBH) of the contacted MS₂.

For MS₂ in contact with graphene at an interlayer rotation angle of zero degrees, the vertical Schottky barrier is absent, which indicates that the heterojunction forms an ohmic contact. The Schottky barriers are n-type with Φ_V = 0.035 eV, 0.021 eV, and 0.042 eV for PtS₂/G at angles of 7.3°, 13.9° and 30°, respectively, and with Φ_V = 0.017 eV, 0.015 eV and 0.057 eV for PtS₂/G at angles of 7.3°, 13.9° and 30°, respectively. Since the values of the n-type SBH are all less than 0.06 eV in all cases, the MS₂/G heterojunction can be regarded as a robust quasi-ohmic contact against the interlayer rotation angles.

3.3 Charge redistribution

To illustrate the nature of charge transference and Fermi level shift at the MS₂-graphene interface, the charge density differences of non-twisted PdS₂/G and PtS₂/G are plotted in Fig. 4(a) and (b). Due to the difference in electronegativity between the sulfur and carbon atoms, an electrostatic potential difference between the two constituents will be established upon the adhesion of graphene on MS₂, leading to a spontaneous polarization at the interface. The generated vertical electric field points from MS₂ to graphene. Such charge transference from PtS₂ (PdS₂) to the graphene layer has been verified by further Bader charge analysis (Fig. 4(c)). Bader charge analysis can specifically describe the influence of the interlayer angle on the interlayer charge transference with a measurable value. The Y-axis Q is defined as the average electronic quantity gained by a single carbon atom, which is shown in the equation below:
where n and N represent the order and total number of carbon atoms of the graphene monolayer in the MS₂/graphene bilayer; Q_G and Q_C represent the electric quantity of a single carbon atom of the calculated system and the electric quantity of a pristine carbon atom, respectively.

Fig. 4 Plot of charge density differences of PdS₂/G (a) and PtS₂/G (b). The green and orange areas represent electron increase and decrease, respectively. The purple, cyan, yellow and gray spheres represent the Pd, Pt, S and C atoms, respectively. (c) Charge variation per C atom versus rotation angles in the MS₂/G systems.

In Fig. 4(c), the calculated Bader charge of the graphene layer is negative, which means that graphene loses electrons and MS₂ is n-type doped. In our work, the relative difference of the Bader charge can quantitatively represent the charge redistribution around the interface. The results show that this charge transference is getting weaker as the rotation angle reaches 30 degrees. As the rotation angle increases, the number of electrons that the graphene layer gains becomes slightly less in both the PtS₂ and PdS₂ cases. This slight change in the number of transferred electrons might be related to the different local interfacial geometries in the heterostructures with different rotation angles.

3.4 Tunneling barrier

The tunneling barrier is an important figure of merit to evaluate the charge injection efficiency of an electrical contact. The tunneling barrier height Φ_TB is defined as the potential energy above E_f between the non-twisted MS₂ and graphene surfaces indicated by the red line in Fig. 5. The tunneling barrier heights Φ_TB of the heterostructures with different rotation angles are listed in Table 1. The value of Φ_TB elevates as the rotation angle increases in the PtS₂/G systems, while this characteristic is completely opposite in the PdS₂/G systems. The tunneling probabilities T_B from graphene to MS₂ are measured using a square potential barrier model such as:
where m is the effective mass of a free electron, ω_B is the width of the tunneling barrier (defined as the full width at half maximum of Φ_TB), and ℏ is the reduced Planck constant. The value of T_B is computed by this equation to be 13.0 (10.0), 16.5 (9.7), 16.5 (13.9), and 22.1 (11.4) % as the rotation angle increases from 0° to 7.3°, 13.9° and 30° for PdS₂/G (PtS₂/G), respectively. These tunneling probabilities are not high compared with the nearly 100% value in some bulk metal contacted systems such as Ni contacted phosphorene⁵⁵ and Ti contacted MoS₂.⁵⁶ Interfacial strong covalent bonds are usually formed in these electrical contacts with high tunneling probabilities. The relatively low tunneling probability in the PtS₂/G and PdS₂/G heterostructures comes from the weak interlayer interaction.

Fig. 5 Plot of electric potential versus z potential for PdS₂/G (a) and PtS₂/G (b). (c) Tunneling barrier of two types of heterojunctions as a function of rotation angle.

4. Conclusions

In conclusion, the structure and electronic properties of PtS₂/G and PdS₂/G heterojunctions with different interlayer rotation angles are studied by DFT calculations. The structure and electronic properties of MS₂/G heterojunctions are almost robust under stacking. The two layers are combined by weak van der Waals interactions, and the heterojunction becomes more stable as the AC direction of graphene coincides with the ZZ direction of MS₂. A slight charge transference from graphene to PtS₂ (PdS₂) occurs, resulting in a vertical build-in internal electric field and PtS₂ (PdS₂) being n-doped. Schottky barriers change in a small range from 0 eV to 0.057 eV in both heterobilayers with different orientations between the two monolayers. In terms of the robust quasi-ohmic nature against the interlayer rotation angle, MS₂/G heterobilayers are potential electrical contacts for high performance nano-electronic devices.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by the National Natural Science Foundation (No. 61675032, 61671085, and 11604019), the National Basic Research Program of China (973 Program) under Grant No. 2014CB643900, the Open Program of the State Key Laboratory of Functional Materials for Informatics, the National Natural Science Foundation for Theoretical Physics special fund “cooperation program” (No. 11547039), and the Shaanxi Institute of Scientific Research Plan projects (No. SLGKYQD2-05).

References

1. K. S. Novoselov, A. K. Geim and S. V. Morozov, Science, 2004, 306, 666–669.
2. A. K. Geim and K. S. Novoselov, Nat. Mater., 2009, 6, 183–191.
3. S. Zhang, J. Zhou and Q. Wang, Proc. Natl. Acad. Sci. U. S. A., 2015, 112(8), 2372–2377.
4. F. Q. Wang, J. Yu and Q. Wang, Carbon, 2016, 105, 424–429.
5. P. Miro, M. Audiffred and T. Heine, Chem. Soc. Rev., 2014, 43, 6537–6554.
6. T. Heine, Acc. Chem. Res., 2015, 48(1), 65–72.
7. H. I. Karunadasa, E. Montalvo and Y. Sun, Science, 2012, 335, 698–702.
8. C. Zhu, X. Mu and P. A. van Aken, Angew. Chem., Int. Ed., 2014, 53, 2152–2156.
9. E. Gourmelon, O. Lignier and H. Hadouda, Sol. Energy Mater. Sol. Cells, 1997, 46, 115–121.
10. B. Hinnemann, P. G. Moses and J. Bonde, J. Am. Chem. Soc., 2005, 127(15), 5308–5309.
11. G. Zhang, H. J. Liu, J. H. Qua and J. H. Li, Energy Environ. Sci., 2012, 9(4), 1190–1209.
12. K. F. Mak, C. Lee, J. Hone, J. Shan and T. F. Heinz, Phys. Rev. Lett., 2010, 105, 136805.
13. B. Radisavljevic, A. Radenovic, J. Brivio, V. Giacometti and A. Kis, Nat. Nanotechnol., 2011, 6(3), 147–150.
14. F. Grønvold and E. Rost, Acta Chem. Scand., 1956, 10, 1620–1634.
15. A. Kjekshus and F. Grønvold, Acta Chem. Scand., 1959, 13, 1767.
16. F. Grønvold, H. Haraldsen, A. Kjekshus and R. Söderquist, Acta Chem. Scand., 1960, 14, 1879–1893.
17. P. Miro, M. Ghorbani-Asl and T. Heine, Angew. Chem., Int. Ed., 2014, 53, 3015–3018.
18. Y. D. Zhao, J. S. Qiao, P. Yu, Z. X. Hu, Z. Y. Lin, S. P. Lau, Z. Liu, W. Ji and Y. Chai, Adv. Mater., 2016, 28, 2399–2407.
19. W. X. Zhang, Z. S. Huang, W. L. Zhang and Y. R. Li, Nano Res., 2014, 7(12), 1731–1737.
20. K. S. Novoselov, A. Mishchenko, A. Carvalho and A. H. Castro Neto, Science, 2016, 353(6298), 1–25.
21. Y. Cai, C. P. Chu, C. M. Wei and M. Y. Chou, Phys. Rev. B: Condens. Matter Mater. Phys., 2013, 88(24), 245408.
22. B. Liu, J. A. Baimova, C. D. Reddy, S. V. Dmitriev, W. K. Law, X. Q. Feng and K. Zhou, Carbon, 2014, 79, 236–244.
23. W. Yang, G. Chen and Z. Shi, et al., Nat. Mater., 2013, 12(9), 792.
24. K. Roy, M. Padmanabhan, S. Goswami, T. P. Sai, G. Ramalingam, S. Raghavan and A. Ghosh, Nat. Nanotechnol., 2013, 8(11), 826–830.
25. W. Zhang, C. P. Chuu and J. K. Huang, et al., Sci. Rep., 2014, 4(7484), 3826.
26. J. E. Padilha, A. Fazzio and A. J. da Silva, Phys. Rev. Lett., 2015, 114(6), 066803.
27. G. Giovannetti, P. A. Khomyakov, G. Brocks, P. J. Kelly and J. van den Brink, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 073103.
28. L. Britnell and L. A. Ponomarenko, Science, 2011, 335, 947.
29. G. H. Lee, Y. J. Yu and X. Cui, et al., ACS Nano, 2013, 7(9), 7931–7936.
30. S. Sattar and U. Schwingenschlögl, ACS Appl. Mater. Interfaces, 2017, 9(18), 15809–15813.
31. J. Pu, Y. Yomogida, K. K. Liu, L. J. Li, Y. Iwasa and T. Takenobu, Nano Lett., 2012, 2(8), 4013–4017.
32. W. J. Yu, Z. Li, H. L. Zhou, Y. Chen, W. Yang, Y. Huang and X. F. Duan, Nat. Mater., 2013, 12(3), 246–252.
33. H. Tian, Z. Tan and C. Wu, et al., Sci. Rep., 2014, 4, 5951.
34. Y. Murata, E. Starodub, B. B. Kappes, C. V. Ciobanu, N. C. Bartelt, K. F. McCarty and S. Kodambaka, Appl. Phys. Lett., 2010, 97(14), 14–16.
35. B. B. Kappes, A. Ebnonnasir, S. Kodambaka and C. V. Ciobanu, Appl. Phys. Lett., 2013, 102(5), 5–10.
36. Y. Murata, S. Nie, A. Ebnonnasir, E. Starodub, B. B. Kappes, K. F. McCarty, C. V. Ciobanu and S. Kodambaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85(20), 205443.
37. Z. Ni, L. Liu, Y. Wang, Z. Zheng, L. J. Li, T. Yu and Z. Shen, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 125404.
38. J. M. B. Lopes dos Santos, N. M. R. Peres and A. H. Castro Neto, Phys. Rev. Lett., 2007, 99, 256802.
39. Y. Wang, Z. Ni, L. Liu, Y. Liu, C. Cong, T. Yu, X. Wang, D. Shen and Z. Shen, ACS Nano, 2010, 4, 4074–4080.
40. J. M. B. Lopes dos Santos, N. M. R. Peres and A. H. Castro Neto, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 86, 155449.
41. Z. Wang, Q. Chen and J. Wang, J. Phys. Chem. C, 2015, 119(9), 4752–4758.
42. W. Kohn and L. J. Sham, Phys. Rev., 1965, 140(4A), A1133–A1138.
43. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6(1), 15–50.
44. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50(24), 17953.
45. G. Kresse and G. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59(3), 1758.
46. D. M. Ceperley and B. J. Alder, Phys. Rev. Lett., 1980, 45(7), 566.
47. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77(18), 3865.
48. S. Grimme, J. Comput. Chem., 2006, 27(15), 1787–1799.
49. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13(12), 5188.
50. A. Kolmogorov and V. Crespi, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 235415.
51. B. Heldt, Nanoscale, 2011, 3(9), 3883.
52. F. Zhang, W. Li, Y. Ma, Y. Tang and X. Dai, RSC Adv., 2017, 7(47), 29350–29356.
53. L. Yu, Y. H. Lee and X. Ling, et al., Nano Lett., 2014, 14(6), 3055.
54. W. Zhang, J. K. Huang and C. H. Chen, et al., Adv. Mater., 2013, 25(25), 3456.
55. Y. Pan, Y. Wang and M. Ye, et al., Chem. Mater., 2016, 28(7), 2100–2109.
56. H. Zhong, R. Quhe and Y. Wang, et al., Sci. Rep., 2016, 6, 21786.
57. J. Gusakov, B. K. Tay and V. Gusakov, Phys. Status Solidi A, 2016, 213(11), 2834–2837.

---