First-principles analysis of seven novel phases of phosphorene with chirality

Abstract

In this work, seven novel phases of phosphorene were predicted to be existent by first-principles calculations, including six kinds of enantiomers corresponding to three kinds of structures with chirality. It is the first time to introduce chirality into two-dimensional (2D) phosphorus. The Poisson’s ratios have been investigated and show normal behavior, rather than the negative one of monolayer or bulk black phosphorus, due to the structures being non-puckered. Phase transformations of these enantiomers have been studied, revealing that there exists the possibility of transformations between them because of the energy barriers being low, which opens doors to possible applications in shape memory devices. This work may inspire new ideas of developing novel applications based on 2D phosphorus nanomaterials.

---

1. Introduction

Nanomaterials have been investigated widely since graphene¹ was discovered. Hence carbon-based structures,^2,3 as the pioneering nanoscale materials, have been abundantly reported. After that, research developments were made about two-dimensional boron nitride,^4,5 an isoelectronic counterpart of graphene. Progress was also made in non-carbon classes of 2D elemental materials such as silicene,⁶ germanene,⁷ stanene,⁸ phosphorene⁹ and arsenene.¹⁰ Transition metal dichalcogenides^11,12 have been a focus for researchers in recent years, which make the 2D family more diversified. In the popular area of 2D nanomaterials, phosphorenes possess many interesting properties such as a negative Poisson ratio,^13,14 strain sensitivity,¹⁵ potential photocatalytic,¹⁶ gas sensing¹⁷ properties and anisotropic intrinsic lattice thermal conductivity.¹⁸ In this work, chirality is introduced into 2D phosphorus for the first time, which may inspire more work to develop this unique property for potential applications.

Most recently, many attempts to develop new structures of phosphorene have been undertaken. Monolayered blue phosphorus (β-P),¹⁹ γ-P and δ-P (ref. 20 and 21) were predicted by Zhu et al. using density functional theory (DFT). Polymorphs of phosphorene with square (ε-P, ζ-P) and pentagonal (η-P, θ-P) units²² were proved to be stable by BOMD (Born–Oppenheimer molecular dynamics) simulation and computation of lattice vibrational spectra. Red phosphorene²³ was designed by analogy with black phosphorene and blue phosphorene. These works made phosphorus more fascinating for nanomaterials science.

Similar methods have been adopted to design new structures from other elements, such as porous boron nitride (PBN),²⁴ which was designed as an analogue of porous graphene,²⁵ and inorganic graphenylene (IGP),^26,27 which was verified to be existent as an analogue of biphenylene carbon (BPC).^28,29 It is mainly because 2D nanostructures are similar in their geometrical conformation that the method of analogy is feasible in predicting new phases. Hence it is possible to design new structures of blue phosphorene by means of analogy with binary elemental 2D structures such as hexagonal BN.

In a previous report²⁶ which used a molecular dynamic simulation, by removing the H atoms of PBN, an analogue of porous graphene, it will spontaneously transform to a new structure named inorganic graphenylene (IGP), a structure similar to BPC. Interestingly, it is feasible to find new structures in 2D phosphorene by replacing the corresponding atoms in PBN and IGP. In this work, the two kinds of atoms in PBN and IGP are substituted by an upper P atom and a lower P atom, respectively. Vibrational spectra and BOMD confirmed their stability. In addition, differently from absolutely one-atom-layered BN, chirality exists in two-atom-layered phosphorus-based 2D materials although they are only a monolayer.

2. Calculation method and models

In this work, simulations were based on density functional theory (DFT). The Vienna ab initio simulation package (VASP)³⁰ was implemented to optimize the structures and investigate their properties. The ion-electron interactions were depicted by projector augmented waves (PAW)³¹ when the function of Perdew, Burke and Ernzerhof (PBE)³² based on the generalized gradient approximation (GGA) was adopted to describe the exchange and correlation potential. We set the cut-off energy to 400 eV in all calculations. To build a 2D structure, a vacuum slab of approximately 20 Å was adopted. A 7 × 7 × 1 Γ-centered Monkhorst-Pack³³ k-point grid was used for sampling of the Brillouin zone for structure optimization and a 15 × 15 × 1 Γ-centered Monkhorst-Pack k-point grid was adopted for band structure calculation. The threshold of convergence was set to 1 × 10⁻⁴ eV and 0.01 eV Å⁻¹ for the self-consistent field (SCF) and ion steps, respectively. All the BOMD simulations were carried out using an NVT ensemble with a 1.0 fs time step and a 5 ps duration.

3. Results and discussion

3.1 Structure and stability of porous phosphorene, IGP-P, DHPP and DHPP-2

There are three kinds of chemical bonds in porous phosphorene, including P–P bonds within hexagonal cycles, P–P bonds connecting neighboring hexagonal cycles and P–H bonds. The optimized geometry structure and band structure of porous phosphorene are shown in Fig. 1(b). The lattice parameters of porous phosphorene are 9.83 Å for a and b, 90° for α and β and 120° for γ. The indirect band gap of porous phosphorene is 2.30 eV, with the lowest point of the CBM (conduction band minimum) located at G and highest point of the VBM (valence band maximum) located at M. To confirm the correctness of the calculations, the lattice parameters and band structure of monolayered blue phosphorene have been calculated (see Fig. 1(a)). It is 3.28 Å and 1.94 eV for the modulation of the lattice vectors and the band gap in our calculation, respectively, which agrees with previous work.¹⁹

Fig. 1 Geometric configurations and band structures of blue phosphorene and four new phases. The arrows indicate the highest point of the VBM and lowest point of the CBM. The upper position and lower position of atoms are denoted by deep colors and light colors, respectively. Unit cells are denoted with a black outline.

To find other new structures, we substituted these two kinds of atoms in IGP with upper P atoms and lower P atoms. This new phase of phosphorene was named inorganic graphenylene based on phosphorus (IGP-P), whose structure and band structure is shown in Fig. 1(c). There are three kinds of bonds in IGP, including P–P bonds in 4- and 6-membered rings, respectively, and in both 4- and 6-membered rings. The lattice parameters of IGP-P are 8.91 Å for a and b, 90° for α and β and 120° for Γ point. The indirect band gap of porous phosphorene is 1.89 eV, with the lowest point of CBM and the highest point of the VBM located at k point.

It is interesting to wonder whether porous phosphorene will spontaneously transform to another phase when the H atoms are removed. Therefore, a BOMD simulation was adopted and resulted in a new structure, which we have named dehydro-porous phosphorene (DHPP), rather than IGP-P. A movie of the BOMD can be found in the ESI.† There are 3-, 4- and 5-membered rings in a unit cell of DHPP. The reason that porous phosphorene doesn’t transform to IGP-P when H atoms are removed is that it is more favorable in energy to stay in the DHPP phase. The cohesive energy is given in Table 1. The cohesive energy represents the energy that would be required to decompose a structure into isolated atoms. A more negative cohesive energy indicates a more stable structure. The formula to calculate the cohesive energy is E_c = E_sys − ∑n_xE_x, where E_c, E_sys, n_x and E_x represent cohesive energy, the energy of a system, the number of x atoms and the energy of one x atom, respectively. Hence it could be seen that all the structures may be stable with respect to thermodynamical stability. The blue phosphorene shows a cohesive energy of −3.47 eV; there is no significant difference compared with those of the new structures in this work. It should be noted that there is no sense in comparing blue phosphorene with porous phosphorene because they are made of different elemental atoms, while comparing blue phosphorene with IGP-P, DHPP and DHPP-2 makes sense. As for an application, the transformation between porous phosphorene and DHPP by removing hydrogen atoms indicates that it may be a good idea to develop porous phosphorene as a hydrogen storage material.³⁴

Table 1 Lattice parameters, cohesive energy, cohesive energy (considering zero point energy) and band gap of blue phosphorene and four new phases

Structure	a (Å)		b (Å)	E c (eV)	E c,z (eV)	E g (eV)
Blue phosphorene		3.28		−3.47	−3.42	1.94
Porous phosphorene		9.83		−2.69	−2.58	2.30
IGP-P		8.91		−3.33	−3.28	1.89
DHPP		9.45		−3.34	−3.29	1.55
DHPP-2	10.42		7.81	−3.33	−3.28	1.70

---

The lattice parameters of DHPP are 9.45 Å for a and b, 90° for α and β and 120° for γ. There is a direct band gap of 1.55 eV at the gamma point, which will be favorable for developing it as a photocatalyst and as a material used for photoelectric conversion without involving energy consumption by phonons. Furthermore, we built another structure similar to DHPP and named it DHPP-2. It has a different lattice from the other structures. The lattice parameters are 10.41 Å and 7.81 Å for a and b, 90° for α and β and 120° for γ. It has an indirect band gap of 1.70 eV, with the highest point of the CBM located between G and K and the lowest point of the VBM located at G.

For more detailed information about the electronic properties of these structures, partial densities of states (PDOS) were calculated as plotted in Fig. 2. For porous phosphorene, H atoms contribute at the CBM but rarely at the VBM, while the density of states from the P atoms contributes to both the VBM and the CBM. For a different quantum number l, p orbitals contribute the most at shallow levels while s orbitals dominate at deep levels lower than −6 eV (vs. the Fermi level). This phenomenon, which appears in all four structures, is in accordance with general quantum mechanics. Hence it is the p orbitals of P atoms that mostly contribute at the VBM and CBM in IGP-P, DHPP and DHPP-2. Band gaps read from the PDOS are consistent with those read from the band structures, although a Gaussian algorithm was used for smearing which made the curves smooth.

Fig. 2 Partial density of states of porous phosphorene (a), IGP-P (b), DHPP (c) and DHPP-2 (d), in which contributions from a different quantum number l and from different kinds of atoms are separated.

To study the vibrational spectra of these phosphorus allotropes, we used the PHONOPY package³⁵ to analyse the phonon dispersion. The phonon dispersion spectra of porous phosphorene, IGP-P, DHPP and DHPP-2 are shown in Fig. 3. There are no significant imaginary frequencies of these three phases, which indicates that they may be stable and possible to synthesize experimentally. Since the electronic energy based on the DFT calculation is related to 0 K, and the entropy correction should be included in order to relate to real conditions, free energy results were used to predict the thermodynamic stability further. Helmholtz free energy was plotted as Fig. 3(e). The free energy F(T) as a function of temperature of blue phosphorene agrees with a previous report.³⁶ Consideration of the zero point energy gives rise to a slightly higher ground state energy for these phases at 0 K. Moreover, it can be seen that while blue phosphorene (which is as stable as black phosphorene with respect to free energy³⁶) is more stable than all the other phases, the F(T) of DHPP is lower than those of IGP-P and DHPP-2 over the whole range of temperatures. It should be noted that comparing the F(T) of porous phosphorene with those of the other phases makes no sense because it includes different kinds of atom, while the other phases only contain phosphorus atoms. Meanwhile, BOMD simulations at 700 K were performed and showed that these four structures were still intact after 5 ps and that the whole frames were vibrating regularly. The movies of these BOMD simulations can be seen in the ESI.†

Fig. 3 Phonon dispersion spectra of (a) porous phosphorene, (b) IGP-P, (c) DHPP and (d) DHPP-2. (e) Helmholtz free energy F(T) as a function of temperature.

3.2 Chirality of IGP-P, DHPP and DHPP-2

Similarly to binary elemental 2D materials, the projection of blue phosphorene in the z direction looks like a honeycomb. However, phosphorene is not an absolutely one-atom-layered structure. There are two atom layers in monolayered phosphorene, so it is possible to involve chirality into systems of phosphorus-based 2D structures. Our calculations concluded that chirality exists in IGP-P, DHPP and DHPP-2 and that each phase has two kinds of enantiomer. The structures of the enantiomers can be seen in Fig. 4. We named these enantiomers with R and S as a suffix. It is easy to identify an R or S enantiomer when the 4-membered ring is examined. If the two upper atoms in a 4-membered ring are at the upper left and bottom right from projection in the z direction, it is an S enantiomer. If the lower atoms in a 4-membered ring are at the upper left and bottom right from projection in the z direction, it is an R enantiomer. In fact, no matter how an enantiomer is rotated, this method to identify them is correct.

Fig. 4 Enantiomers of IGP-P, DHPP and DHPP-2. (a) IGP-P-S, (b) IGP-P-R, (c) DHPP-S, (d) DHPP-R, (e) DHPP-2-S and (f) DHPP-2-R. The upper position and lower position of atoms are denoted by deep colors and light colors, respectively. (g) Differences of adsorption energy corresponding to enantiomers of CHFClBr, alanine and glutamic acid adsorbed on IGP-P, in which the solid circles relate to the left axis and the hollow circles relate to the right axis.

When chirality is considered, new applications such as chiral catalysis, molecular self-assembly and polarized optics will be feasible in the future. Many previous works were conducted to utilize the chirality of nanomaterials for various applications.^37,38 As for the chemical properties, it is significant that chirality plays an important role. To adsorb and separate different enantiomers of molecules, it could be possible to prepare these novel nanomaterials as a functional resin. By modifying the functional groups, based on their intrinsic chirality, chiral catalysis may be achieved for reactions of biochemical molecules including drugs and imaging agents, considering the good compatibility of phosphorus in biochemistry. However, as for the physical properties, currently the cohesive energy and band gap of the two kinds of enantiomers are the same. Therefore, in this work, all the calculations are based on the S enantiomers.

To emphasize the importance of their chirality, the adsorption of chiral molecules on IGP-P-S was adopted as an example. There is only one position of P atoms in IGP-P-S, and around a P atom there in turn exist 12-, 4- and 6-membered rings clockwise. Hence each neighboring P atom of this P is different after we picked it. When enantiomers of molecules with chirality are adsorbed on this P atom, the interaction of the substrate and the three different ligands of the molecule is different. It is mostly van der Waals forces that contribute the difference, and a DFT-D3 correction^39,40 was adopted. Three kinds of molecules containing different number of atoms including CHFClBr, alanine and glutamic were tested. It is well known that the bigger a molecule is, the stronger the van der Waals force effects will be. For the smallest molecule, which only has five atoms, CHFClBr shows the smallest difference between its enantiomers (see Fig. 4g). In contrast, glutamic acid shows the biggest effect of chiral adsorption. It follows that this kind of chiral material may show selective adsorption of chiral molecules, which would be useful for applications. Detailed data can be seen in the ESI.†

3.3 Poisson ratio of porous phosphorene, IGP-P and DHPP

To study the Poisson ratio of porous BP, IGP-P and DHPP, we calculated the strain along different directions at the same time, a method which was adopted in early work.^13,14 For calculating the Poisson ratio, a square lattice was adopted. The transformational matrix of the lattice vectors is presented as . The vectors matrix after redefinition is obtained as , where a₂, b₂ and c₂ are for vectors after redefinition and a₁, b₁ and c₁ are for vectors before redefinition. The redefined lattices are shown in Fig. 5(c) and (f), corresponding to IGP-P and DHPP.

Fig. 5 Surface plot of (ε_x, ε_y) and the corresponding strain energies of (a) IGP-P and (d) DHPP. Surface plot of (ε_x, ε_y) and the corresponding ε_z of (b) IGP-P and (e) DHPP. Redefined lattice adopted for the Poisson ratio’s calculation of (c) IGP-P and (f) DHPP.

We set strain along the x and y directions in a range of −0.02 to 0.02 and with an interval of 0.004; data points can be seen in Fig. 5 as black dots. The strain energy is defined as E_s = E_x,y − E_0,0, where E_s, E_x,y and E_0,0 represent the strain energy, the total energy of the system when strain is adopted and the total energy of the system when no strain is adopted in two directions. The relationship of E_s, ε_x and ε_y is fitted in the following polynomial equation: E_s = aε_x² + bε_y² + cε_xε_y, in which a, b and c are fitting parameters. The values of a, b and c can be seen in Table 2. Stress along the x(y) direction was calculated by the formulation σ_x(y) = V₀⁻¹∂E_s/∂ε_x(y), where V₀ is a constant as the equilibrium volume. Uniaxial deformations along the x and y directions are indicated by black lines shown in Fig. 5, which are denoted by σ_x = 0 and σ_y = 0, respectively. The corresponding Poisson’s ratios are evaluated as ν_y = c/2a and ν_x = c/2b, where ν_y and ν_x are defined as ν_y = −ε_x/ε_y and ν_x = −ε_y/ε_x, respectively. The same method was adopted in an earlier report¹⁴ and we successfully repeated the results of it.

Table 2 Fitting numbers and calculated Poisson ratios of IGP-P and DHPP

Structure	a	b	c	ν y	ν x
IGP-P	174.22	194.61	132.21	0.379	0.340
DHPP	159.73	122.40	50.89	0.159	0.208

---

The Poisson ratios calculated are 0.379 (0.159) and 0.340 (0.208) for ν_y and ν_x for IGP-P (DHPP), respectively. It seems that the Poisson ratios of IGP-P are larger than those of DHPP in both in-plane directions. This phenomenon is mainly due to the cell volume of DHPP being larger than that of IGP-P, which makes the force of the atoms in the DHPP lattice less than those in IGP-P when uniaxial strain is adopted. Hence there is less influence on lattice deformation in another direction.

The strain along the z direction was calculated and is plotted in Fig. 5(b) and (e) for IGP-P and DHPP, respectively. Previous work¹⁴ shows that both bulk and monolayered black phosphorus exhibit negative Poisson’s ratios along the armchair direction in response to perpendicular uniaxial strains. However, in this work it is not the same. Under uniaxial strain along the x(y) direction, strain along the z direction is increased (reduced) as the sheet is compressed (stretched), in which a positive strain indicates stretching and negative strain indicates compressing. It is different from the existence of a negative Poisson’s ratio in black phosphorene. It is the structure that the properties depend on. To be specific, the puckered structure of black phosphorene doesn’t exist in blue phosphorene, where the puckered structure means that when you look at the structure from the top of black phosphorene, the lower atoms would be blocked by upper atoms. However, this phenomenon doesn’t appear in the blue one. Hence it follows that IGP-P and DHPP exhibit normal Poisson ratio behavior in response to uniaxial deformations oriented in all directions.

3.4 Phase transformation of enantiomers of IGP-P and DHPP

The phase transformations of different kinds of structures were investigated in previous research.^19,20,22 It is interesting as to whether these new phases in this work are possible to transform to each other or not. For DHPP and IGP-P, their cell shapes are the same while their cell volumes are different. Actually, the cell volume could be presented as lattice constant |a| for a hexagonal lattice. We have plotted the energy of each phase when lattice constant |a| varies in Fig. 6(a). It seems that IGP-P could be more stable than DHPP at a relatively smaller cell volume, which could be approached by increasing the pressure. Hence transformation from DHPP to IGP-P under higher pressure is feasible from their energy dependence on cell volume.

Fig. 6 (a) Transformation from IGP-P to DHPP under external pressure, which is revealed to be feasible by their energy dependence on cell volume, denoted by lattice constant |a|. (b) Transformation path and energy barrier from DHPP-S to IGP-P-S. (c) Transformation path and energy barrier between two enantiomers of DHPP. (d) Whole scale of transformation of enantiomers of IGP-P and DHPP, where E_a represents the activation energy of transformation.

The transition paths of enantiomers of IGP-P and DHPP were also studied in this work. A generalized solid-state nudged elastic band (GSS-NEB) method^41–43 was adopted to find the energy barriers of their transformation. It is unlikely that the IGP-P-S transforms to IGP-P-R directly, while it is easy for enantiomers of DHPP to be transformed between each other. If an enantiomer of IGP-P must transform to another, the transition path is as follows. At first, a transformation to DHPP with the same chirality would take place, then a transformation to DHPP with different chirality. Finally, a transformation to IGP-P with the different chirality would occur. The paths of the transformations are shown in Fig. 6(b) for DHPP-S to IGP-P-S and (c) for DHPP-S to DHPP-R. The path is bilaterally symmetric for the plotted transformation path of different enantiomers of the same phase, while it is bilaterally asymmetric for the plotted transformation path of different phases with the same chirality. The calculated potential barriers of IGP-P-S to DHPP-S (IGP-P-R to DHPP-R) and DHPP-S to DHPP-R (DHPP-R to DHPP-S) are 0.066 eV per atom and 0.077 eV per atom, respectively. Such a small energy barrier indicates that it is possible to obtain racemic IGP-P and DHPP when the system absorbs enough energy from the environment. It is reasonable to expect that if some pure IGP-S was obtained, it would transform to a mixture of IGP-P-S, IGP-P-R, DHPP-S and DHPP-R by suitable heating.

It is likely that the GSS-NEB method underestimates the energy barrier, so it may be a little higher in reality. We also tested another method used in previous research,^19,22 in which lattice constants of the initial state were redefined to coincide with the lattice of the final state. The calculations show that there is no significant difference between these two methods due to the similar lattice configuration.

4. Conclusion and prospects

In this work, four new phases of phosphorene were predicted to exist; three of these have the property of chirality, and enantiomers. The transformations of them have been investigated and it seems to be possible to obtain racemic production in certain conditions. Their Poisson’s ratios are found to be normal in the z direction, rather than the negative one of black phosphorene, due to non-puckered structures. Their property of chirality could possibly be used in optical materials and chiral catalysis. Transformation between porous phosphorene and DHPP by removing hydrogen atoms indicates that it may be a good application to develop porous phosphorene as a hydrogen storage material.

For experimental fabrication,⁴⁴ top-down methods and bottom-up methods are both possible ways to obtain these large-hole structures. When black phosphorus is cleaved by large-volume ions in solution environments, the ions are embedded into bulk layers. Hence the layers are separated by the force of the ions. In this procedure, defects of large holes may be generated. As for bottom-up methods, it depends on the substrates when chemical vapor deposition (CVD) is adopted. If specific substrates were used, these materials could be synthesized.

In addition, the catalytic activity of phosphorene has not been thoroughly studied so far. The Lewis base property makes phosphorene likely to be a good catalyst in some specified reactions. Involving the chirality, chiral catalysis by phosphorus-based nanomaterials will be interesting. More experimental and theoretical study may stimulate the potential of these novel materials.

Acknowledgements

This work was supported by the State Key Program of Natural Science of Tianjin (Grant No. 13JCZDJC26800), the MOE Innovation Team (IRT13022) of China, the National Natural Science Foundation of China (Grants No. 142100, 21433008, 91545106), and the foundation of State Key Laboratory of Coal Conversion (Grant No. J15-16-908).

References

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva and A. A. Firsov, Science, 2004, 306, 666–669.
2. H. Zhang, W. Ding and D. K. Aidun, J. Nanosci. Nanotechnol., 2015, 15, 1660–1668.
3. H. Zhang, W. Ding, K. He and M. Li, Nanoscale Res. Lett., 2010, 5, 1264–1271.
4. D. Pacile, J. C. Meyer, C. O. Girit and A. Zettl, Appl. Phys. Lett., 2008, 92, 133107.
5. X.-F. Zhou, X. Dong, A. R. Oganov, Q. Zhu, Y. Tian and H.-T. Wang, Phys. Rev. Lett., 2014, 112, 085502.
6. P. De Padova, C. Quaresima, C. Ottaviani, P. M. Sheverdyaeva, P. Moras, C. Carbone, D. Topwal, B. Olivieri, A. Kara, H. Oughaddou, B. Aufray and G. Le Lay, Appl. Phys. Lett., 2010, 96, 261905.
7. Z. Ni, Q. Liu, K. Tang, J. Zheng, J. Zhou, R. Qin, Z. Gao, D. Yu and J. Lu, Nano Lett., 2012, 12, 113–118.
8. P. Tang, P. Chen, W. Cao, H. Huang, S. Cahangirov, L. Xian, Y. Xu, S.-C. Zhang, W. Duan and A. Rubio, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 90, 121408.
9. H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tomanek and P. D. Ye, ACS Nano, 2014, 8, 4033–4041.
10. S. Zhang, Z. Yan, Y. Li, Z. Chen and H. Zeng, Angew. Chem., Int. Ed., 2015, 54, 3112–3115.
11. P. K. Chow, E. Singh, B. C. Viana, J. Gao, J. Luo, J. Li, Z. Lin, A. L. Elias, Y. Shi, Z. Wang, M. Terrones and N. Koratkar, ACS Nano, 2015, 9, 3023–3031.
12. Y. Ding and B. Xiao, RSC Adv., 2015, 5, 18391–18400.
13. J.-W. Jiang and H. S. Park, Nat. Commun., 2014, 5, 4727.
14. M. Elahi, K. Khaliji, S. M. Tabatabaei, M. Pourfath and R. Asgari, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 91, 115412.
15. R. Fei and L. Yang, Nano Lett., 2014, 14, 2884–2889.
16. B. Sa, Y.-L. Li, J. Qi, R. Ahuja and Z. Sun, J. Phys. Chem. C, 2014, 118, 26560–26568.
17. L. Kou, T. Frauenheim and C. Chen, J. Phys. Chem. Lett., 2014, 5, 2675–2681.
18. R. Fei, A. Faghaninia, R. Soklaski, J.-A. Yan, C. Lo and L. Yang, Nano Lett., 2014, 14, 6393–6399.
19. Z. Zhu and D. Tomanek, Phys. Rev. Lett., 2014, 112, 176802.
20. J. Guan, Z. Zhu and D. Tomanek, Phys. Rev. Lett., 2014, 113, 046804.
21. J. Guan, Z. Zhu and D. Tomanek, ACS Nano, 2014, 8, 12763–12768.
22. M. Wu, H. Fu, L. Zhou, K. Yao and X. C. Zeng, Nano Lett., 2015, 15, 3557–3562.
23. T. Zhao, C. Y. He, S. Y. Ma, K. W. Zhang, X. Y. Peng, G. F. Xie and J. X. Zhong, J. Phys.: Condens. Matter, 2015, 27, 265301.
24. Y. Ding, Y. Wang, S. Shi and W. Tang, J. Phys. Chem. C, 2011, 115, 5334–5343.
25. P. Kuhn, A. Forget, D. Su, A. Thomas and M. Antonietti, J. Am. Chem. Soc., 2008, 130, 13333–13337.
26. E. Perim, R. Paupitz, P. A. S. Autreto and D. S. Galvao, J. Phys. Chem. C, 2014, 118, 23670–23674.
27. J. R. Feng and G. C. Wang, Comput. Mater. Sci., 2016, 111, 366–373.
28. G. Brunetto, P. A. S. Autreto, L. D. Machado, B. I. Santos, R. P. B. dos Santos and D. S. Galvao, J. Phys. Chem. C, 2012, 116, 12810–12813.
29. Y. Li, D. Datta, Z. Li and V. B. Shenoy, Comput. Mater. Sci., 2014, 83, 212–216.
30. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
31. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
32. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
33. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
34. J. Luo, H.-B. Zhou, Y.-L. Liu, L.-J. Gui, S. Jin, Y. Zhang and G.-H. Lu, J. Phys.: Condens. Matter, 2011, 23, 135501.
35. A. Togo, F. Oba and I. Tanaka, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 134106.
36. Y. Aierken, D. Cakir, C. Sevik and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 081408.
37. J. L. Zhang, T. C. Niu, A. T. S. Wee and W. Chen, Phys. Chem. Chem. Phys., 2013, 15, 12414–12427.
38. T. Chen, Q. Chen, X. Zhang, D. Wang and L.-J. Wan, J. Am. Chem. Soc., 2010, 132, 5598–5599.
39. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
40. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
41. G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000, 113, 9901–9904.
42. G. Henkelman and H. Jonsson, J. Chem. Phys., 2000, 113, 9978–9985.
43. D. Sheppard, P. Xiao, W. Chemelewski, D. D. Johnson and G. Henkelman, J. Chem. Phys., 2012, 136, 034103.
44. L. Kou, C. Chen and S. C. Smith, J. Phys. Chem. Lett., 2015, 6, 2794–2805.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra25361f

---