Engineering the electronic bandgaps and band edge positions in carbon-substituted 2D boron nitride: a first-principles investigation

Abstract

Modification of graphene to open a robust gap in its electronic spectrum is essential for its use in field effect transistors and photochemistry applications. Inspired by recent experimental success in the preparation of homogeneous alloys of graphene and boron nitride (BN), we consider here engineering the electronic structure and bandgap of C_2xB_1−xN_1−x alloys via both compositional and configurational modification. We start from the BN end-member, which already has a large bandgap, and then show that (a) the bandgap can in principle be reduced to about 2 eV with moderate substitution of C (x < 0.25); and (b) the electronic structure of C_2xB_1−xN_1−x can be further tuned not only with composition x, but also with the configuration adopted by C substituents in the BN matrix. Our analysis, based on accurate screened hybrid functional calculations, provides a clear understanding of the correlation found between the bandgap and the level of aggregation of C atoms: the bandgap decreases most when the C atoms are maximally isolated, and increases with aggregation of C atoms due to the formation of bonding and anti-bonding bands associated with hybridization of occupied and empty defect states. We determine the location of valence and conduction band edges relative to vacuum and discuss the implications on the potential use of 2D C_2xB_1−xN_1−x alloys in photocatalytic applications. Finally, we assess the thermodynamic limitations on the formation of these alloys using a cluster expansion model derived from first-principles.

---

1 Introduction

Hexagonal boron nitride (h-BN) has a bulk crystal structure analogous to that of graphite. An isolated layer of h-BN, with a two-dimensional (2D) honeycomb structure, is thus analogous to graphene.^1,2 But despite the structural equivalence, the electronic properties of these 2D materials differ greatly: graphene (referred to as G here) is a zero-bandgap semiconductor, whereas a h-BN monolayer (referred to as BN here) is an insulator with a wide band gap E_g > 5 eV. The simultaneous similarity in crystal structure (with a relatively small lattice mismatch of ∼1.8%) and contrast in electronic behavior, offer potential for a number of exciting applications based on combinations of the two materials. One current direction of experimental research in this area is the controlled synthesis of in-plane junctions between conductive G and insulating BN, with the purpose of developing one-atom-thick integrated circuits,^3–5 where it is obviously necessary to prevent the unintentional intermixing of the two phases. On the other hand, one might want to intentionally mix G and BN in order to achieve intermediate electronic bandgaps that could be useful in electronic or optical devices. However, there are significant thermodynamic constraints for the formation of homogeneous G–BN alloys as a strong driving force exists to segregate G and BN domains/nanophases.⁶ Some theoretical work has therefore focused on understanding the effect of domain distribution on the electronic structure^7–12 and mechanical properties^13,14 of the mixed system. In a recent article, Lu et al. demonstrated the synthesis of highly homogeneous G–BN alloys supported on ruthenium.¹⁵ These authors found that the energetics of mixing and demixing processes are modified by the presence of the metal support. These alloys are still metastable with respect to phase separation, but their synthesis can be achieved under non-equilibrium conditions at high temperatures, followed by rapid quenching to prevent the diffusion of species towards segregated domains. This experimental progress calls for a better theoretical understanding of the properties of highly homogeneous G–BN alloys, where the formation of domains is inhibited.

In the present work, we study how the electronic structure of 2D carbon/boron nitride alloys is determined not only by composition, but also by the distribution of the ions at a given composition. This is in the spirit of previous research which has revealed interesting and potentially tuneable variations of semiconductor bandgaps within the configurational space of ion distributions at a fixed composition.^16–19 We start from the BN end-member of the solid solution, which already has a large gap, and study the behavior of the gap upon C substitution. In addition to the bandgap, we consider the alignment of the valence and conduction band edges with respect to the vacuum reference level; this alignment is important in understanding the electronic behaviour of the interfaces that these alloys form with other metallic or semiconducting materials, and also the potential of these nanostructures for photochemistry applications. The compositions considered here can be written as C_2xB_1−xN_1−x, i.e., C replaces the same number of B and N atoms, keeping the ratio B/N = 1, and low values of x (x < 1/3). These compositions give bandgap values in a range that is useful for optoelectronic and photochemistry applications.

2 Methods

For the electronic structure calculations, the 2D materials were simulated using periodic slabs, with layers separated by a vacuum gap of fixed width (20 Å) along the c direction. In the lateral directions, the simulation supercell consisted of 3 × 3 unit cells and contained 18 sites (9 B and 9 N sites for pure BN). The cell was substituted with 2, 4 and 6 C atoms to simulate C_2xB_1−xN_1−x compositions with x = 0.11, 0.22 and 0.33, respectively. The symmetrically inequivalent configurations at each composition in the 3 × 3 supercell were generated using the Site Occupancy Disorder (SOD) program.^20,21 The criterion for two configurations to be symmetrically equivalent is that an isometric transformation converts one into the other, and the list of possible transformations is obtained from the space group of the parent structure, in combination with the supercell translations.

Our first-principles calculations were based on the density functional theory (DFT) as implemented in the Vienna Ab initio Simulation Package (VASP).^22–25 The projector augmented wave (PAW) method^26,27 was used to describe the interaction between ionic cores (including the 1s level on each atom) and valence electrons. An energy cutoff of 520 eV was used for the plane wave basis set expansion. Integrations in the k-space were made using a 8 × 8 × 1 uniform mesh of points within the reciprocal lattice of the supercell (which corresponds to a 24 × 24 × 1 mesh in the reciprocal lattice of the unit cell). The atomic positions and lateral lattice parameters (the c parameter was kept constant) were optimized with the Perdew–Burke–Ernzerhof (PBE)²⁸ exchange correlation functional, which is based on the generalized gradient approximated (GGA). At the final geometries, single-point calculations based on the screened hybrid functional of Heyd, Scuseria and Ernzerhof (HSE06)^29,30 were performed to obtain the electronic structure. The HSE06 functional is known to yield accurate predictions of electronic bandgaps in semiconductors,³¹ in contrast to the typical underestimation resulting from GGA functionals.

The evaluation of configuration energies in a larger supercell for the thermodynamic analysis were performed using a cluster expansion model,³² including both nearest and next nearest neighbour clusters for pairs and triplets. The interaction parameters were fitted to DFT energies in the smaller cell.

3 Results and discussion

3.1 Bandgaps

The supercell with composition C₂B₈N₈ (x = 0.11) has only 3 symmetrically different site-occupancy configurations: one configuration where the two C atoms are in nearest-neighbor sites, forming a C–C bond or ‘dimer’, and two configurations without C–C bonds. The supercell with composition C₄B₇N₇ (x = 0.22) has 30 different configurations, with varying degrees of C aggregation, ranging from the cases where all four C atoms form a cluster (tetramer), to cases where each C atom is “isolated” (no C–C bonds). Our first observation is that configurations with higher C content tend to have narrower bandgaps, as expected from the electronic structures of pure G and BN. However, there are also significant variations in bandgaps amongst configurations corresponding to the same compositions. We have plotted the bandgap versus the average size of C clusters in each configuration (refer to Fig. 1). Clearly, structures with higher carbon aggregation tend to have wider band gaps as compared to those with smaller cluster sizes or isolated C atoms. In the C₄B₇N₇ configuration with the narrowest bandgap (≈2 eV) all the four C atoms are isolated.

Fig. 1 Variation of bandgaps with average cluster sizes in the C₂B₈N₈ cell composition (x = 0.11; black circles), and in the C₄B₇N₇ cell composition (x = 0.22; red squares). The bandgap decreases with C concentration and increases with the average cluster size in the configuration. The bandgap of h-BN (x = 0; green diamond) is given for reference.

We note that a range of bandgap values still exist for a given average cluster size, but this dispersion can also be explained by the C distribution. For example, in the C₄B₇N₇ (x = 0.22) composition, when the average cluster size is 2, there are two possible types of configurations according to the cluster size distribution: configurations have either one triplet and one isolated carbon (“1+3”), or two dimers of carbon atoms (“2+2”). We find that the 1+3 configurations give bandgaps smaller than those of the 2+2 configurations.

In order to understand the origin of the trend described above, we have studied the contribution of the orbitals from different atomic species to the electronic density of states for two C₂B₈N₈ (x = 0.11) configurations: one where the C atoms are isolated, and one where they are forming a C–C dimer (see Fig. 2). We have aligned the electron energies in both cases with respect to the vacuum level, which was determined from the electron potential in the middle of the vacuum gap in the periodic simulation cell. The double substitution (C at N and C at B) gives rise to two defect states, one above the valence band maximum of pure BN and and one below the conduction band minimum of pure BN, in agreement with a previous report.⁸

Fig. 2 Electronic density of states and charge density plot of Highest Occupied Molecular Orbital (HOMO) and Lowest Unoccupied Molecular Orbital (LUMO) for configurations with isolated (“1+1”) C atoms and dimers (“2”) in C₂B₈N₈ (x = 0.22). (a) Density of states, and charge density plots of (b) HOMO and (c) LUMO states of “1+1” configuration. (d) Density of states, and charge density plots of (e) HOMO and (f) LUMO states of “2” configuration. (g) Schematic of the formation of bonding and antibonding states from the C/B and C/N defect states in “1+1” configuration. E_F denotes the Fermi energy and the dotted line at 0 is the vacuum level. Here, C = yellow, B = green and N = blue.

When the C substituents are separate from each other (Fig. 2(a)), the electron density corresponding to the highest occupied band in the substituted system is clearly localized at the C atom occupying the N site (Fig. 2(b)), whereas the density corresponding to the first empty band is localized at the C atom occupying the B site (Fig. 2(c)). However, when the two C atoms form a bond, the defect states (Fig. 2(d)) arise from the mixing of the orbitals from the two C atoms. The highest occupied state in this case is the “bonding” combination (Fig. 2(e)), and the lowest empty state is the “anti-bonding” combination (Fig. 2(f)). Because the energy difference between the “bonding” and “antibonding” states is larger than that between the isolated C/N and C/B defect states, the bandgap increases when dimerization occurs (as shown schematically in Fig. 2(g)). The formation of a bonding defect state also explains why the dimerized configuration is strongly stabilized with respect to the configurations with only C monomers: the calculated energy difference is 1.66 eV. We will return to the discussion of stabilities below.

3.2 Band alignment

The alignment of the band structure with reference to the vacuum level is important as a basis to understand the electronic behavior of interfaces involving these 2D materials. Such referencing also allows us to explore the potential activity in photocatalytic applications, for example, for hydrogen production via water splitting. The possibility of making these nanostructures active for photocatalytic hydrogen production is very attractive because, if achieved, it could lead to simultaneous hydrogen production and storage in combined C/BN systems.

Assuming a single-semiconductor cell configuration, and that other materials properties are satisfied, an ideal photocatalyst for water splitting would have certain important characteristics in its electronic structure. On the one hand, the positions of the conduction and valence band edges should straddle the redox potentials for water photolysis.^33–35 This means that the conduction band edge should be above the energy corresponding to the hydrogen evolution reaction (HER), and the valence band edge should be below the energy of the oxygen evolution reaction (OER), implying also that the bandgap must be wider than 1.23 eV (difference between the HER and the OER). In fact, once loss mechanisms are accounted for, a bandgap of 2 eV or more is generally necessary.³⁴ On the other hand, the bandgap should not be too wide, in order to allow the adsorption of photons from the visible part of solar radiation. It is known that, in the vacuum scale and at pH = 0, the HER level is located at −4.44 eV, and the OER level is located at −5.67 eV.³⁶ At temperature T and pH > 0, these energy levels are shifted up by (k_BT ln 10) × pH, where k_B is Boltzmann's constant. By referencing the electronic levels in our semiconductor solid solutions to the vacuum level (taken here as the electron potential at the middle of the vacuum gap), we can assess whether the band edges of the semiconductor are in a favorable position to catalyze the solar splitting of water under a given set of conditions.

The positions of the band edges of C₂B₈N₈ (x = 0.11) and C₄B₇N_0.7 (x = 0.22) configurations with respect to the vacuum potential and to the HER and OER electrode potentials at pH = 0 and pH = 7 are shown in Fig. 3(a) and (b) respectively. Configurations are sorted in each case by their total energies (from lower to higher), and the band edge positions of pure BN are given for reference. Pure BN has has band edges straddling the HER and OER redox potentials (configurations with this property will be referred to as straddling configurations), and therefore may exhibit photocatalytic activity, but its bandgap is too wide for efficient solar energy absorption. For the C₂B₈N₈ supercell composition there is one configuration with the two C atoms forming a dimer (configuration 1), and two configurations (2 and 3) with isolated C substitutions (no C–C bonds). The dimerized configuration has band edges straddling the HER and OER redox potentials at both pH = 0 and pH = 7, but its bandgap (around 4.1 eV) is still too wide for efficient use of solar energy. The configurations with isolated C atoms have narrower bandgaps of around 3.3 eV, but their valence band minimum is too high, above the OER level at pH = 0, and therefore would not be able to donate a hole for the OER reaction to occur. At pH = 7 the valence band edges for the isolated C configurations fall just below the OER level, but are still too near it to drive this reaction efficiently.

Fig. 3 Bandgaps and band edge positions calculated with the HSE06 functional for the symmetrically inequivalent configurations of C_2xB_1−xN_1−x, for (a) C₂B₈N₈ (x = 0.11), (b) C₄B₇N₇ (x = 0.22), and (c) for configurations with only dimers in C₆B₆N₆ (x = 0.33) with respect to the vacuum potential. The configurations are arranged on the x axis in increasing order of their total energies. The energy levels for the hydrogen evolution reaction (HER) and oxygen evolution reaction (OER) are represented by blue and green lines, respectively, both at pH = 0 (solid line) and at pH = 7 (dashed line).

For the C₄B₇N₇ composition (Fig. 3(b)), the 30 symmetrically different configurations can be divided into four groups, in order of stability. In the first group (configurations 1 and 2 in Fig. 3(b)), with the lowest total energies, the configurations consist of maximally segregated carbon distributions, i.e. forming C₄ clusters. The bandgap for both these configurations is approximately 3.3 eV, and their valence band maximum falls slightly above the OER level at pH = 0, and slightly below it at pH = 7. Configurations following these (second group, configurations 3 to 9 in Fig. 3(b)) are either of “2+2” type, i.e., containing two isolated C–C dimers per supercell, or of “1+3” type (mixed, i.e. one triplet and one isolated C). In particular, for configurations 7 and 9, which are both of “2+2” type, the valence band maximum falls below the OER potential and therefore the conduction and valence band edges straddle the redox potentials for water photolysis (even at pH = 0). Their bandgap (around 3.7 eV) is reduced with respect to the dimerized configurations in the C₂B₈N₈ composition, but is still too wide for efficient photocatalysis. In the third group (configurations 10 to 26 in Fig. 3(b)), the configurations are mixtures of “1+1+2” and “1+3”, and exhibit narrower band gaps. Lastly (fourth part, 27–30), the configurations consist entirely of isolated C substitutions (“1+1+1+1”), and have band gaps in the range of 2.2 to 2.5 eV, which are suitable for efficient absorption of solar light. But again, the isolated C configurations have a valence band maximum above the OER level and therefore would not be suitable for a single-semiconductor photolysis cell not even at higher pH values (pH = 7). However, it is interesting to note that semiconductors with such band positions can in principle be utilized as a cathode, in conjunction with other materials (such as TiO₂) for the anode, in a heterojunction photocatalyst. To be useful as a photocathode, the semiconductor band edges are only required to straddle the HER level.^33–35

Finally, we explore whether we can still achieve configurations straddling both the OER and the HER level by increasing the C concentration (assuming that such high concentrations can be achieved in a homogeneous solid solution, which will be discussed in the next section). We have seen before that only configurations with dimerized substitutions have a valence band edge low enough to allow hole donation for the OER reaction. Therefore, we have considered all the symmetrically distinct “2+2+2” configurations with the C₆B₆N₆ (x = 0.33) composition. Our results show that not all these are straddling configurations at pH = 0 (Fig. 3(c)). However, at pH = 7, all configurations straddle the OER and HER levels, although the most stable configuration (number 1) has the valence band maximum too close to the OER level. With the bands positions in the correct region of the energy scale, it is possible to achieve the desired band alignment if the reaction conditions are changed to slightly higher pH or temperatures.

Therefore, our band alignment results show that (a) the band structure of these alloys could be engineered to make them work as single-semiconductor photocatalyst for the water splitting reaction under conditions near to room temperature and neutral pH, if the C substituents are distributed forming C–C dimers, but these configurations tend to have too wide bandgaps for efficient solar energy utilization; and (b) configurations with only isolated C substituents, if successfully synthesized at a reasonably high C concentration, could achieve narrower bandgaps in the desired range for solar light absortpion but with a valence band edge too high for single-semiconductor water splitting photocatalysts. These configurations would satisfy the band edge position requirements for a cathode in a heterojunction photocatalyst.

3.3 Stability

We briefly discuss here the relative energies of the configurations, and also the stability of the alloy with respect to phase separation. The calculated mixing energies (energies of the alloy with respect to the corresponding amounts of pure BN and graphene) are positive and high for all the compositions and configurations, confirming that the alloy formation is always a highly endothermic process (Fig. 4). For reference, we have also plotted the average energy of the alloy as a function of composition (seen as red line in Fig. 4), as calculated in a larger (6 × 6) cell for each composition. This was done using a cluster expansion³² model with parameters fitted to DFT energies in the smaller cell; the agreement between the energies calculated at DFT level and those obtained from the cluster expansion is illustrated in the inset of Fig. 4. The configurations evaluated with the cluster expansion and included in the average energy all have the C_2xB_1−xN_1−x stoichiometry.

Fig. 4 Crosses show the DFT mixing energies (per formula unit) for C_2xB_1−xN_1−x configurations and the red line represents the average mixing energies calculated for a converged sample of random configurations using a cluster expansion model. Inset shows the mixing energies predicted by the model versus DFT for C_2xB_1−xN_1−x.

The very high mixing energies shown in Fig. 4 indicate that the free energy of mixing, despite the stabilizing effect of the configurational entropy, will also be positive at any significant C concentration. To illustrate this, consider the maximum configurational entropy for this type of alloy, using the ideal expresssion:

S = −2k_B(x ln x + (1−x)ln(1−x)). (1)

where the factor of 2 comes from the two sublattices with mixed occupancy. At x = 0.1 and T = 1000 K, −ST = −56 meV per formula unit, which is more than 5 times smaller in absolute value than the mixing energy at that composition.

We explained above how the formation of C–C bonds lowers the total energy with respect to two isolated C atoms; we see here that the same trend is observed for larger clusters, which are always energetically favoured over smaller ones. Thermodynamics clearly drives the system towards the segregation of BN and C domains, and against the formation of an alloy, as has been discussed previously by other authors.^6,9,15 Our results show that isolated carbon substitutions are always the least stable, and are above the average energy line, particularly for higher concentrations, which suggest that such configurations might be difficult to achieve. It is important to note that these alloys must always be synthetized by a non-equilibrium procedure^15,37 and therefore it is difficult to predict the extent to which the energy differences between configurations are relevant. However, our results suggest that in a metastable sample with significant C substitution, obtained by quenching from high temperature, isolated C atoms would occur rarely, if at all. A sample with only isolated C substitutions could only be obtained by such methods at very low C concentrations, where the configurational entropy dominates the distribution. At higher C concentration, a more sophisticated synthesis technique allowing atom-by-atom control would be required. On the other hand, configurations with isolated C–C dimers, which are also potentially interesting in photocatalysis, might be easier to synthesize because the C–C interaction stabilizes them with respect to isolated C substitutions. However, they are still metastable with respect to phase segregation.

4 Conclusions

Our work elucidates the relationship between the distribution of C substituents in 2D boron nitride and the resulting electronic structure. We show that the bandgap decreases with increasing C concentration, as expected, and that the reduction in gap is more pronounced the more disperse the distribution of the substituents is. Hence the smallest bandgaps for a given composition correspond to configurations with only isolated C substituents, that is, with no C–C bonds. For example, we have found that the substitution of isolated C atoms at a concentration of only x = 0.22 brings the bandgap down to 2 eV.

From the point of view of potential applications in electronics, it is important to note that this small gap arises from narrow bands associated with defect states, leading to large effective masses and low mobilities. However, for larger concentrations of C substitutions the edge bands become wider. Therefore when engineering the bandgap of this system one would need to find a balance between gap opening and mobility.

The analysis of the electronic structure shows that the clustering of the C substituents leads to mixing of the wavefunctions corresponding to these defect states, forming bonding and anti-bonding levels. This effect is responsible for the increase in the band gap with C clustering.

We have also reported the variation of the band edge positions with respect to composition and ion distribution. We found that configurations consisting of C–C dimers have band alignments which are favourable for single-semiconductor photocatalysis of the water splitting reaction under conditions near to room temperature and neutral pH, but these configurations tend to have too wide bandgaps for efficient solar energy utilization. On the other hand, configurations with only isolated C atoms have smaller bandgaps. Their band edge positions do not straddle the OER level, but satisfy the requirements for a cathode in a heterojunction photocatalyst (straddling the HER level). However, our thermodynamic analysis shows that configurations with isolated C atoms would be very difficult to obtain due to the strong tendency of C substituents to cluster in the BN matrix.

Acknowledgements

The authors thank the UK-India Education and Research Initiative (UKIERI) for funding this collaborative project through a Trilateral Partnership Grant IND/CONT/2013-14/054. Via the UK's HPC Materials Chemistry Consortium, which is funded by EPSRC (EP/L000202), this work made use of ARCHER, the UK's national high-performance computing service.

References

1. J. C. Meyer, A. Chuvilin, G. Algara-Siller, J. Biskupek and U. Kaiser, Nano Lett., 2009, 9, 2683–2689.
2. A. Nag, K. Raidongia, K. P. S. S. Hembram, R. Datta, U. V. Waghmare and C. N. R. Rao, ACS Nano, 2010, 4, 1539–1544.
3. P. Sutter, R. Cortes, J. Lahiri and E. Sutter, Nano Lett., 2012, 12, 4869–4874.
4. M. P. Levendorf, C. Kim, L. Brown, P. Y. Huang, R. W. Havener, D. A. Muller and J. Park, Nature, 2012, 488, 627–632.
5. Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Lei, X. Yang, J. Zhang, J. Yu and K. P. Hackenberg, et al. , Nat. Nanotechnol., 2013, 8, 119–124.
6. L. Ci, L. Song, C. Jin, D. Jariwala, D. Wu, Y. Li, A. Srivastava, Z. Wang, K. Storr and L. Balicas, et al. , Nat. Mater., 2010, 9, 430–435.
7. J. da Rocha Martins and H. Chacham, ACS Nano, 2011, 5, 385–393.
8. M. Bernardi, M. Palummo and J. C. Grossman, Phys. Rev. Lett., 2012, 108, 226805.
9. J. Zhu, S. Bhandary, B. Sanyal and H. Ottosson, J. Phys. Chem. C, 2011, 115, 10264–10271.
10. A. K. Manna and S. K. Pati, J. Phys. Chem. C, 2011, 115, 10842–10850.
11. N. Kumar, K. Moses, K. Pramoda, S. N. Shirodkar, A. K. Mishra, U. V. Waghmare, A. Sundaresan and C. N. R. Rao, J. Mater. Chem. A, 2013, 1, 5806–5821.
12. Q. Peng and S. De, Phys. E, 2012, 44, 1662–1666.
13. Q. Peng, A. R. Zamiri, W. Ji and S. De, Acta Mech., 2012, 223, 2591–2596.
14. Q. Peng, X.-J. Chen, W. Ji and S. De, Adv. Eng. Mater., 2013, 15, 718–727.
15. J. Lu, K. Zhang, X. F. Liu, H. Zhang, T. C. Sum, A. H. C. Neto and K. P. Loh, Nat. Commun., 2013, 4, 2681.
16. P. Piquini, P. A. Graf and A. Zunger, Phys. Rev. Lett., 2008, 100, 186403.
17. Y. Seminovski, P. Palacios, P. Wahnón and R. Grau-Crespo, Appl. Phys. Lett., 2012, 100, 102112.
18. D. O. Scanlon and A. Walsh, Appl. Phys. Lett., 2012, 100, 251911.
19. D. Santos-Carballal, A. Roldan, R. Grau-Crespo and N. H. de Leeuw, Phys. Rev. B: Condens. Matter Mater. Phys., 2015 , in press, http://arxiv.org/abs/1504.00268.
20. R. Grau-Crespo, S. Hamad, C. R. A. Catlow and N. H. de Leeuw, J. Phys.: Condens. Matter, 2007, 19, 256201.
21. R. Grau-Crespo and U. V. Waghmare, in Molecular Modeling for the Design of Novel Performance Chemicals and Materials, ed. B. Rai, CRC Press, 2012, ch. 11, pp. 303–326.
22. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
23. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
24. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
25. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
26. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
27. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
28. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
29. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
30. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2006, 124, 219906.
31. J. Heyd, J. E. Peralta, G. E. Scuseria and R. L. Martin, J. Chem. Phys., 2005, 123, 174101.
32. C. Wolverton and D. de Fontaine, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 8627–8642.
33. H. Wang, L. Zhang, Z. Chen, J. Hu, S. Li, Z. Wang, J. Liu and X. Wang, Chem. Soc. Rev., 2014, 43, 5234–5244.
34. M. G. Walter, E. L. Warren, J. R. McKone, S. W. Boettcher, Q. Mi, E. A. Santori and N. S. Lewis, Chem. Rev., 2010, 110, 6446–6473.
35. T. Hisatomi, J. Kubota and K. Domen, Chem. Soc. Rev., 2014, 43, 7520–7535.
36. S. Trasatti, Pure Appl. Chem., 1986, 58, 955–966.
37. K. Raidongia, A. Nag, K. P. S. S. Hembram, U. V. Waghmare, R. Datta and C. N. R. Rao, Chem. – Eur. J., 2010, 16, 149–157.

---

Footnote

† Present address: School of Engineering and Applied Sciences, Harvard University, Cambridge, Massachusetts 02138, USA.

---