Exploring the stability and electronic structure of beryllium and sulphur co-doped graphene: a first principles study

Abstract

First principles density functional theory (DFT) calculations have been performed to explore the stability, structural and electronic properties of Be and S co-doped graphene sheets. The band-gap of graphene has been tuned by co-doping with beryllium and sulphur at different sites. The results show that by co-doping graphene with Be and S, the band-gap increases from zero up to 0.58 eV depending on the doping sites. The cohesive and formation energies of the systems were also determined. All the isomers formed by exploring different doping sites differ notably in stability, bond length and band-gap. Nevertheless, the planar structure of all the systems investigated was preserved even after geometry optimisation. The majority of the isomers that correspond to co-doping at non-equivalent sites favour higher band-gap opening, but lower stability, than the other set of isomers with equivalent doping sites. Bader charge analysis was adopted to account for charge distribution in the systems. As a result of the difference in electronegativity among carbon atoms and the impurities, it was observed that electron accumulation occurred more on the carbon atoms in the proximity of Be and S than at any other position in the graphitic systems investigated.

---

Introduction

Graphene is actively investigated as a possible replacement for silicon for use in micro, and even nanoelectronics. In silicon-based electronics, useful devices are made by doping with electropositive and electronegative elements. Similarly, doping will be needed to create useful devices in graphene-based electronics. This article uses DFT to investigate the atomic and electronic structure of one type of doping of graphene, namely substitution of carbon atoms by atoms of beryllium and sulphur.

Graphene is a two dimensional sp² hybridized carbon atom material with a honey comb lattice and the unit cell has two carbon atoms in non-equivalent positions which form two interpenetrating triangular sublattices (see Fig. 1). It is known to be a fundamental and integral part of other forms of graphitic carbon such as fullerenes (bucky ball), carbon nanotubes and graphite.^1–4 While its existence was theoretically predicted by P. R. Wallace many years back,¹ it was not believed at that time to be physically attainable as a 2D crystalline material because it was adjudged to be thermodynamically unstable.⁵ In 2004 however, Novoselov et al.⁶ experimentally isolated the graphene crystal through micromechanical exfoliation of pyrolytic graphite. Ever since, the material has drawn a lot of attention among researchers in the scientific community. More importantly, graphene has been lauded as a potential material in semiconductor nanoelectronics as a result of its stable hybrid structure, excellent electrical and optical properties and high charge carrier mobility.^3,6–8 In fact the charge carrier mobility, which is about 10⁶ cm² V⁻¹ s⁻¹, has been reported to be 2–3 times higher than that of conventional semiconductors.⁹

Fig. 1 Graphene has two carbon atoms in the unit cell. The two atoms are at non-equivalent sites denoted by α and β where a⃑ and b⃑ symbolized the primitive unit vectors.

Pristine graphene however has zero band gap which makes it extremely difficult for direct electronic device applications. For instance in field-effect transistors (FETs) the role of a band gap in the drain material is crucial without which it would be difficult to switch off the device. To overcome this shortcoming, different techniques can be used to fabricate sizable band gaps in graphene sheets. Some of the techniques that could be employed to open up an energy gap in graphene are: superstructures fabrication, surface modification with atoms or polar molecules,^10–14 chemical doping,^15–22 application of electric field,^23,24 and deposition of graphene on epitaxial substrates.^25–27

Superstructure fabrication like quantum-dot,^28,29 graphene nanoribbons^30,31 or nanomeshes³² induce an energy gap in a material through quantum confinement. However, engineering of such structures dimensional-wise still poses a serious challenge. The design of such superstructures with uniform widths and edges below 10 nm is still quite tedious to attain with top-down techniques such as etching and lithography.³³ Although homogenous structures can be achieved easily with bottom-up approaches,³⁴ pulling off patterned ribbons still remains a contending procedure.

A facile approach to induce an energy gap in graphene is by depositing it on epitaxial substrates like SiC,³⁵ Al₂O₃,³⁶ etc. Apart from the difficulties involved in controlling the morphology and the surface energies of the material, the band gap induced by this technique is not tuneable.³⁷ However, for a tuneable band gap inducing technique, application of electric field could be a viable option.

The application of an external electric field (E) perpendicular to the plane of AB stacked bilayer graphene can induce a band gap tuneable up to 0.25 eV through the breaking of the material inversion symmetry. The carrier mobility is not substantially affected in the process. Despite the robustness of this technique for AB stacked bilayer graphene, the electric field does not favour a band-gap opening directly in monolayer graphene.³⁷ This is simply because the application of this field to monolayer graphene does not affect the inherent symmetry of the material sublattices.

Furthermore, unlike the electric field scheme which is applicable exclusively to few-layer graphene, chemical doping method via substitution of foreign atomic elements (heteroatoms) in place of the carbon atoms in graphene structure has been reported to be a suitable technique for creating sizeable and varying band gaps in monolayer graphene.³⁸ It involves using elements (like N and B) with comparable atomic radius (but different number of valence electrons) as carbon to replace the carbon atoms of graphene. The technique of substitution came into prominence shortly after the isolation of a single graphene layer in 2004, and ever since the method has widely been cited in different works in literature.³⁹

In this study, we focused primarily on chemical doping as an avenue of modifying the electronic structure of graphene for band gap engineering, and the following are the highlights of the previous studies objectified by this technique. For examples, using ab initio calculations, Lee et al.⁴⁰ studied the adsorption of beryllium (Be) on fullerenes and their potential to store hydrogen molecules. They found that the clustering of Be on pristine fullerenes is energetically stable but causes the dissociation of hydrogen molecules. Ferro et al.⁴¹ carried out DFT studies of the absorption, diffusion of Be in graphite and the formation of Be₂C. They reported that, in high quantity, the absorption of Be in graphite can result into formation of Be₂C. Denis et al.⁴² performed first principle studies of S and P-doped graphene. They found that below 0.5% impurity concentration, S-doped graphene has a band gap in the range of 0.1–0.2 eV while P-doped graphene has a higher gap between 0.3 and 0.4 eV. Recently, Ullah et al.³³ investigated structural and electronic properties of Be and Be/B dual doped graphene via density functional theory calculations, and remarked that with a rectangular doping configuration a maximum band gap of 1.44 eV can be obtained for Be-doped graphene at 12.5% impurity concentration. Moreover a maximum band gap of 0.99 eV can be attained for Be–B co-doped graphene at the same concentration. More recently, Huang et al.⁴³ reported electronic structure with optical properties of boron-doped, sulphur-doped and boron and sulphur dual doped graphene. They observed that doping modified the electronic and the optical properties of the doped graphene samples in question, and that the size of the energy gap depends on the impurities concentration in a doped graphene sample. Some other works relating to sulphur-doped graphene or SWCNTs can be found in ref. 44.

It is worth noting that in spite of the numerous previous studies, there has been none that addresses structural properties and electronic structure of Be and S co-doped graphene. Moreover, Be–S co-doping of graphene is expected to make a better semi-conducting material than either Be-doped or S-doped graphene. This is because doping of graphene with sulphur could lead to a material with metallic behaviour,¹⁵ while doping with beryllium could form a material with Fermi-level lying within the valence band,³³ which to us might form a material with lower carrier mobility than Be–S co-doped graphene which has higher Fermi-level but lies within the energy-gap. In view of this, using first principle calculations within the frame work of density functional theory (DFT) we investigated structural properties and electronic structures of Be–S co-doped graphene. We have tried not only to accomplish band gap engineering and tuning in the sheet through Be and sulphur substitutional co-doping, but we have also explored the stability and the dependability of the energy gap on the doping sites across the sublattices in the graphene crystal. In addition, Bader charge analysis^45–47 has been adopted to account for charge distribution in the material under study.

Methods

First-principles numerical non-spin polarised calculations were carried out using ab initio density functional theory as implemented in MedeA VASP software.^48–51 In the calculations, the projected augmented wave (PAW) was used to describe the system of electron–ion interactions, while the generalized gradient approximation (GGA) with PBE prescription⁵² was utilized for exchange correlation. The kinetic energy cut-off for the wave function expansion was set to 400 eV with energy convergence criteria for all the systems considered set to 10⁻⁵ eV. Graphitic systems of 32 atoms with supercell model slab of 4 × 4 were employed in all the calculations. To avoid interlayer interaction between two successive layers of graphene in the simulations, a vacuum spacing of 14 Å was used in between the layers with applied periodic boundary condition along the graphitic plane. All the structures were geometrically relaxed until the Hellmann–Feynman forces converged below 0.02 eV Å⁻¹. Based on the symmetry of the graphitic systems considered and in order to also actualise self-consistent field calculations, a Γ-centred grid of 6 × 6 × 1 k mesh with Gaussian smearing was employed to sample the Brillouin zone. Moreover, a denser k-point grid of 17 × 17 × 1 was adopted when calculating the density of states (DOS). Cohesive energy, E_coh calculations were performed using the following formula:
here x represents C, Be or S; E_x and E_t stand for the atomic energies and electronic ground state energy of the whole system, respectively. The variable N stands for the total number of atoms in the system while n_x is the number of atoms of type x. The negative sign outside the bracket implies that all the stable structures should have positive cohesive energies, otherwise the system is unstable and the magnitude of the value is a measure of the bonding strength.

To calculate the atomic energies of the atoms, we created a 10 Å × 10 Å × 10 Å super cell and centred an atom in the middle of the box. We set the K-integration to Gaussian and performed a single point spin polarised calculation with Γ–k-point. In this case spin polarised calculations was performed because carbon and sulphur atoms are both open shell systems.

Furthermore, we have calculated the difference between the formation energy (FE) of Be–S co-doped and the sum of the formation energies of Be-doped and S-doped graphene to check if this study is in agreement with the recent study by Denis et al.⁵³ who have shown that the FE of certain classes of dual doped graphene is less than the sum of the formation energies of the corresponding mono-doped graphene systems. To calculate the formation energy of co-doped and mono-doped graphene systems, we employed the following formula:

FE_(BeS) = E_BeS + 2μ_C − E_graphene − μ_Be − μ_S

FE_(Be) = E_Be + 2μ_C − E_graphene − μ_Be

FE_(S) = E_S + 2μ_C − E_graphene − μ_S

here E_BeS, E_Be and E_S represent the electronic ground state energy of Be–S co-doped graphene, Be and S doped graphene respectively. E_graphene is the energy of pristine graphene. For the chemical potentials (μ_C, μ_Be and μ_S), the atomic energies of carbon, beryllium and sulphur were used appropriately for the calculation of the formation energy of the systems. To this end, it is worth mentioning that the conclusion drawn from the analysis of the formation energy of dual doped graphene is independent of the sources of the chemical potentials used in the calculation.⁵³

Note: having initially performed spin and non-spin polarised calculations for the geometry optimization of one of the configurations of Be and S co-doped graphene and realised that the energies of the two calculations are equal, along with zero net magnetic moment, we subsequently considered only non-spin polarised calculations for the rest of the configurations in order to save ample amount of time. Besides, Be and S are not the only dopants that preserve non-magnetic nature of graphene. Recently, Hussain et al.⁵⁴ also reported Be and N co-doped graphene to be non-magnetic.

Results and discussion

First, a geometry optimization was carried out on a pure graphene sheet to allow all the carbon atoms to relax. After the optimization, the lattice constant of pristine graphene was found to be 2.46 Å, while the C–C bond length measures as 1.42 Å, as shown in Fig. 2a. These values are in excellent agreement with some reported theoretical³³ and experimental⁵⁵ results. The calculated band structure of the optimized pristine structure of graphene is presented in Fig. 2b which also corroborates previous findings in relations to zero band gap and linear dispersion crossing at the Dirac point.⁵⁶

Fig. 2 (a) The optimised geometry of a 4 × 4 supercell of pristine graphene sheet. (b) The calculated band structure of pristine graphene.

Subsequently, a pristine graphene sheet was co-doped with beryllium and sulphur atoms. Different configurations of the above mentioned impurities in the graphene sheet were considered in order to exploit the effect of doping sites on the electronic properties and stability of the co-doped graphene sheets.

In defining a configuration of a Be–S co-doped graphene (BeSG) system, we substituted two carbon atoms at equivalent sites (i.e. at O and C1) and non-equivalent sites (i.e. at O and C) of pristine graphene for beryllium (blue) and sulphur atom (yellow) as shown in Fig. 3. The co-doped structure with sulphur and beryllium now occupying O and C1 positions respectively is termed the OC1 configuration or simply the OC1 isomer while the OC isomer symbolizes another configuration with S and Be being substituted for two carbon atoms at the sites O and C of the pristine graphene.

Fig. 3 (a) OC1 (which simply refers to as OC1 isomer) exemplifies a configuration with two carbon atoms at same sublattices (i.e. O and C1) substituted for Be and S. The blue colour in the diagram stands for Be while the yellow colour represents sulphur. (b) OC symbolises a configuration with two carbon atoms at different sublattices (i.e. O and C) substituted for Be and S.

Furthermore, in a similar manner to OC and OC1, other isomers of Be and S co-doped graphene (BeSG) system were also generated such that the relative positions of Be and S have been changed successively. Fig. 4 shows an overview of other configurations considered while exploring structural and electronic properties of Be and S co-doped graphene. Configurations OA1, OB1, OC1, OD1, etc. correspond to same sublattices or equivalent sites substitution of Be and S co-doping whereas, OA, OB, OC, or OD matches co-doping of monolayer graphene at different sublattices (or non-equivalent sites). It is worth noting that all the configurations were chosen with respect to the position of S atom at the lattice point, O in the systems.

Fig. 4 Different isomers of Be–S co-doped graphene sheets. The relative position of Be has been updated through positions A, A1, B, B1, C, C1, etc. The sites A1, B1, C1, etc. denote Be atom position at the same sublattices with respect to S atom while other positions of Be atoms are non-equivalent with respect to S-atom. The blue colour represents Be-atom whereas yellow stands for sulphur atom.

Different isomers of Be–S co-doped graphene considered in this study yielded different electronic band structures and density of states. In the following discussion, for the sake of conciseness, we have considered in detail the cohesive energy and formation energy; charge transfer; structural and electronic properties of OD, OB1 and OB isomers of Be–S co-doped graphene sheets.

OD isomer of BeSG

Two carbon atoms of graphene were replaced with Be and S-atom as shown in Fig. 5a. This was followed by a geometry optimisation of the system until it reaches the level of the required accuracy. After the optimisation, it was observed that the planar structure of the system was preserved except for the palpable changes in the adjoining bond lengths. The adjoining bond lengths of S–C and Be–C expand to 1.62 Å and 1.55 Å respectively; this is as a result of the sizes of molecular covalent radii of Be (0.90 Å) and S (1.02 Å) which are larger than the covalent radius of the carbon atom (0.77 Å). Consequently, the C–C bond length nearest to sulphur reduces to 1.37–1.39 Å while that of beryllium reduces to 1.36–1.38 Å (see Fig. 5a). The sizes of the bond lengths of Be–C and C–C reported are in agreement with previous studies.^33,54 Besides, the cohesive energy which gives a measure of the respective stabilities of Be–S co-doped graphene systems was calculated and found to be equal to 7.32 eV per atom for this system with OD isomer. This value is less than 9.19 eV per atom that we calculated for 4 × 4 pristine graphene. Meanwhile, the cohesive energy that we obtained for pristine graphene is in agreement with the previous reports by Ullah et al.³³ (9.21 eV per atom) and Rani and Jindal³⁸ (9.20 eV per atom). Note we have dropped the sign of the cohesive energy based on the model adopted as explained under Methods.

Fig. 5 (a) Optimized geometry of a Be–S co-doped graphene sheet. (b) Total DOS of BeSG with OD configuration. (c) The band structure of the optimised structure with OD configuration. (d) Projected DOS of carbon in BeSG with OD configuration. (e) Projected DOS of Be in BeSG with OD configuration. (f) Projected DOS of sulphur in BeSG with OD configuration.

Bader charge analysis was invoked to calculate the charge transfer among Be, S and Cs in the system. Due to the differing electronegativity among these three atoms, the valence charge of 1.57e (Be) plus 0.26e (S) is transferred to the carbon atoms in the system; the transferred charge is redistributed among the carbon atoms with the charges on each stretching from 3.81e to 4.47e. Moreover, as expected higher value of the charge is transferred to the carbon atoms in the neighbourhood of the impurities.

The calculations of density of states (DOS) and partial density of states (PDOS) were performed to study the effect of the co-doping on the system. The results are shown in Fig. 5 which shows that p_xy (i.e. 2s, 2p_x and 2p_y orbitals hybridization) orbitals of sulphur and beryllium hybridized sturdily with the p_xy of the carbon atoms of BeSG system. The stout overlapping between p_xy orbitals of carbon atoms and sulphur atom (in the range −8.0 to −8.5 eV in the valence region) partly leads to p–σ bonds formation. Likewise, π-bond in the system could partly be attributed to the very weak overlapping of p_z orbitals of the carbon atoms and sulphur atom in the system. As for the beryllium atom, it contributed little to the DOS in the conduction band region. However, while the p_xy orbitals of Be hybridized strongly with p_xy orbitals of the carbon atoms around the Fermi level between 0 and 1.0 eV, the p_z orbital overlapped weakly with that of the carbon atoms contributing to π-bond formation of the system.

The electronic band structure of the system was also calculated. A direct band-gap of 0.58 eV is induced around the Dirac point (see Fig. 5c). Here the position of the Fermi-level remains as it is in pristine graphene. This is as a result of isoelectronic nature of BeS co-coped graphene with respect to pristine graphene. That is, equal proportion of holes and electrons are introduced into carbon network simultaneously. S is electron rich and has two electrons more than C while Be is electrons deficient and has two electrons less than C as such the net change in the number of electrons in the system is zero. Our results show that the band-gap induced by Be and S co-doped graphene with OD configuration is greater than those reported in the previous studies for single Be³³ or S⁴³ doping of graphene. Interestingly, a transistor with a graphene platform requires a minimum energy gap of 0.4 eV to operate suitably in on or off mode.⁵⁷ Thus the value (0.58 eV) from our calculation meets this specification. In addition, DFT usually underestimates band gaps hence the actual value for the band gap in this calculation might be greater than the value we reported. This significant result spurred us to perform further calculations on BeSG with different configurations.

OB1 isomer of BeSG

In the following case, S and Be are substituted accordingly for two carbons of graphene at O and B1 positions which are equivalent sites (see Fig. 6a). After the substitution, geometry optimisation was performed which still leads to the retention of the planar shape of the system but there were changes in the adjoining bond lengths of the structure. The optimised structure of the new system is presented in Fig. 6a. The bond length of S–C increases (relatively to C–C of pristine graphene) and falls within the range 1.54–1.62 Å while that of Be–C lies within the range 1.53–1.54 Å. Due to the large sizes of the covalent atomic radii of the aforementioned impurities relative to a carbon atom, the C–C bonds in the neighbourhood of S reduced to values within the range 1.36–1.40 Å whereas the values within the range of 1.36–1.38 Å (Fig. 6a) was observed for the bond lengths of C–C in the proximity of the Be impurity in the system under study. Unlike OD isomer, there is fluctuation in the bond lengths of S–C and Be–C of BeSG with OB1 isomer. For this system, the cohesive energy was calculated to be 7.35 eV per atom which is higher than that of OD isomer.

Fig. 6 (a) The optimised structure of OB1 isomer of Be–S doped graphene. (b) The band structure of the optimised structure with OB1 configuration. Total DOS of Be–S with OB1 configuration. (c) The band structure of the optimised structure with OB1 configuration. (d) Projected DOS of carbon in Be–S with OB1 configuration. (e) Projected DOS of Be in Be–S with OB1 configuration. (f) Projected DOS of sulphur in Be–S with OB1 configuration.

Bader charge analysis reveals that the valence charge of 1.57e (Be) along with 0.48e from S is transferred to the carbon atoms in the system. It was noted that the size of the charge transferred from sulphur was greater than the corresponding transfer in the previous system with OD isomer. After the transfer of charges from the impurities to the system, charge reformation on carbon atoms occurs and stretches from 3.80e to 4.61e. In addition, we found as expected that more charges were transferred to the carbon atoms adjacent to the impurities in the graphene sheet.

To study the role of dopants in this system with this aforementioned configuration, the calculations of DOS and PDOS were performed. The results are displayed in Fig. 6. It is evident from the plots that the DOS contribution to the conduction band originated mainly from p_xy orbital of C atoms and little contribution from S. There is also minor contribution from p_z of S to these bands. For the valence bands, besides C, the population of this band is from p_xy orbitals of both S and Be. p_xy of C and Be strongly overlapped and form a σ bond which covers an energy range of −2.0 to 0.0 eV. This wider range of energy of interaction between p_xy of C and Be over the corresponding overlap in OD isomer might account for the slight increase in the cohesive energy of this system. Similarly, σ bond is also formed between p_xy orbitals of C and S. In this case the interaction covers an energy range from −8.5 eV to −3.5 eV, although it is more intense between −8.5 eV to 8.0 eV.

The electronic band structure calculations show a band gap opening of 0.22 eV around the Fermi level of this system (see Fig. 6b). The gap opening around the Fermi level reflects the isoelectronic nature of Be–S co-doping. This value is less than the observed band gap of the previous structure with OD isomer. Thus the difference in the size of the band gaps can be ascribed to the configuration of the dopants in the given structures. This variation in the band gaps of the two isomers of the dual doped graphene due to the relative positions of the substituted dopants in them is consistent with the past studies^38,58 on B and N co-doped graphene.

OB isomer of Be–S co-doped graphene

When Be and S are substituted for two carbon atoms of graphene at two non-equivalent sites (i.e. at O and B) nearest to each other (see Fig. 7a), similar effects to the previously discussed isomers are observed after the optimisation of the structure. The dopants interact with the other carbon atoms through sp² hybridization. After optimisation of the structure, the bond lengths of Be–S, S–C and Be–C were found to be 1.69 Å 1.62 Å and 1.50 Å respectively. These bond lengths are greater than the optimised C–C bond length of graphene. Consequently the immediate C–C bond lengths adjacent to the dopants are reduced leading to the distortion of some of the adjoining hexagonal carbon rings of the system. Though there were significant changes in the adjoining bond lengths of the atoms of the structure, the planar shape of the system is preserved.

Fig. 7 (a) Optimized geometry of a Be–S co-doped graphene sheet. (b) Total DOS of BeSG with OB configuration. (c) The band structure of the optimised structure with OB configuration. (d) Projected DOS of carbon in BeSG with OB configuration. (e) Projected DOS of Be in BeSG with OB configuration. (f) Projected DOS of sulphur in BeSG with OB configuration.

For this isomer the cohesive energy of the system is equal to 7.41 eV per atom which is the highest value among the isomers of Be–S co-doped systems we considered, thus the most stable. Analogously, this configuration shows conformity with the most stable configuration of nitrogen and boron co-doped graphene as cited by Nath et al.⁵⁹ The detailed explanation of this favourable configuration has been provided in the immediate section.

Bader charge analysis confirms that the valence charge of 1.56e from Be along with 0.20e from S is transferred to the carbon atoms in the system. For this isomer, the size of the charge transferred from sulphur is less than the corresponding transfer in the previous system with OB1 isomer. After the transfer of charges from the impurities to the system, charge redistribution on carbon atoms occurs which ranges from 3.75e to 4.52e. Moreover just like the previous cases, we observed that more charges were transferred to the carbon atoms adjacent to the impurities in the graphene sheet.

The results from the DOS and PDOS are presented in Fig. 7. The DOS contribution around the Dirac point comes mainly from p_xy of C, S and less from Be. There is also minor contribution from p_z orbital of S. The p_xy orbitals of Be does not only strongly hybridized with p_xy of C atoms over the range of −2.3 to 1.0 eV but also overlapped relatively softly over an interval of −4 to −9 eV and this leads to p–σ bonds near the Fermi levels. Similarly, S and C which strappingly hybridized over a wider energy range (−14 to −1.0 eV) also forms p–σ bond. The π bond between C and the dopants is very weak. This is illustrated by the slight overlapping between p_z of carbon and that of the impurities. As it can be seen from the band structures in Fig. 7c a gap opening of 0.28 eV which deviates slightly away from the Dirac point is observed. The deviation of the extremum of the gap of this system away from the Dirac point might be attributed to the loss of hexagonal symmetry of some of the carbon rings in the vicinity of the impurities in the system and as such the linear dispersion relationship at the Dirac point is lost. The non-uniformity in the hexagonal symmetry of the system in question is exemplified by the relative larger bond lengths of Be–S and S–C which have been uncovered in the preceding paragraphs under this isomer.

Formation energies of BeS co-doped graphene

We have included this section to show that this work is in line with the recent study by Denis et al.⁵³ where it has been shown that certain classes of dual doped graphene have tendency to be more stable and easier to synthesise than their corresponding mono-doped graphene. For example, they pointed out that Al–X (where X = B, N, O) dual doped graphene are more stable than Al-doped graphene, and as such they would be easier to synthesise. In light of their report, we have calculated the formation energies of Be–S co-doped graphene involving all the isomers in this study and compared it with the formation energies of Be-doped and S-doped graphene. The result is presented in Table 1. We would like to reiterate once more that for the calculation of FE, atomic energies of the atoms were considered for the chemical potentials.

Table 1 Differences between the formation energies (FE) of 4 × 4 BeS co-doped graphene and the sum of FE of 4 × 4 Be and S singly doped graphene^a

Configuration	FE(BeS) (eV)	∑(FE(Be) + FE(S)) (eV)	FE(BS) − ∑(FE(Be) + FE(S)) (eV)
a Note: FE(Be) = 2.09 eV and FE(S) = 3.10 eV.
OA	2.45	5.19	−2.74
OA1	2.39	5.19	−2.80
OB	0.61	5.19	−4.58
OB1	2.39	5.19	−2.80
OC	2.45	5.19	−2.74
OC1	3.80	5.19	−1.39
OD	3.47	5.19	−1.72
OQ	3.47	5.19	−1.72

---

From Table 1, it can be seen that 0.61 eV which is associated with OB isomer has the least formation energy, followed by 2.39 eV that is shared by both the OA1 and the OB1 configuration. On the other hand, OC1 isomer has the highest FE which implies that is the least stable among the configurations considered in this study.

The above values show that Be and S atoms prefer to replace C–C bond of graphene since other configurations have higher FE. These favourable sites for Be–S co-doped graphene have also been echoed in recent studies^54,60 to be suitable sites for certain classes of dual doped graphene. So far, the only exception to this proposition is SiB along with GeB dual doped graphene. It is interesting to point out that the FE of BeS with OB1 isomer is lower than FE of Be or S-doped graphene. Similarly, OA, OA1, OB1 and OC have lower FE than S-doped graphene. Meanwhile, previous reports^61,62 showed S-doped graphene has been successfully synthesised. Therefore we could generalise that it would be easier to synthesise Be–S co-doped graphene with these configuration than to synthesise S-doped graphene. In addition Table 1 reveals that the FE computed for different isomers of Be–S co-doped graphene is less than the FE resulting from the sum of Be and S-doped graphene. This shows that there is synergy between Be and S which lowered the FE significantly than the sum of FE of Be and S singly doped graphene.

Summary

Band gap opening is important for developing electronic devices with graphene platforms. Nevertheless, it is not enough to create a finite band gap in graphene without a method of tuning energy gaps to suit different applications. Chemical doping of graphene at varying concentration via atomic substitution has constantly been reported as a good technique to actualise a band-gap tuning in graphene. However, doping of graphene with atomic impurities at increasing concentration leads to decrease in the cohesive energy of the material, and thus reduced stability. In order to overcome this shortcoming, the technique in this study which involves Be–S co-doping of graphene across the sublattices, is apparently a viable option of achieving band-gap tuning while keeping the impurities concentration in check, and thus the stability of the system. With this background, the plot of the variation of the cohesive energies of Be–S co-doped graphene with the different configurations is shown in Fig. 8a and supported with the Table 2 which includes the summary of all the results of the calculations.

Fig. 8 (a) The variation of the cohesive energies of Be–S co-doped graphene with the different configurations. (b) The variation of the energy gaps of Be–S co-doped graphene with the different configurations.

Table 2 Summary of the configurations against the corresponding cohesive energies (coh), band-gaps and difference between formation energies (FE) of 4 × 4 BeSG and the sum of FE of 4 × 4 Be and S singly doped graphene^a

Configuration	Cohesive energy (eV per atom)	FE(BS) − ∑(FE(Be) + FE(S)) (eV)	Band gap (eV)
a Note: FE(Be) = 2.09 eV and FE(S) = 3.10 eV, coh of pristine graphene = 9.19 eV.
OA	7.35	−2.74	0.45
OA1	7.35	−2.80	0.18
OB	7.41	−4.58	0.28
OB1	7.35	−2.80	0.18
OC	7.35	−2.74	0.20
OC1	7.31	−1.39	0.31
OD	7.32	−1.72	0.58
OQ	7.32	−1.72	0.39

---

It is evident from the plot that the pair of dopants, Be and S, substituted as nearest neighbour (OB isomer) in graphene is the most stable configuration since the cohesive energy is the highest while the formation energy the least among the configurations considered. This specific configuration (i.e. isomer OB) yields 7.41 eV per dopant whereas other isomers have cohesive energies that lie within 7.31–7.35 eV per dopant.

In addition, the band structures have been calculated for co-doped graphene systems via sampling route Γ → M → K → Γ → A → L → H → A → L → M → K → H of the hexagonal Brillouin zone. Certain patterns of variations are observed in the electronic band structures and in the DOS of all the co-doped systems considered. It was found that pristine graphene exhibits a linear band at k-point of its hexagonal BZ with zero band-gap, whereas a finite band-gap has been observed in each of the aforementioned isomers of Be–S co-doped graphene systems investigated. Fig. 8b reveals the energy gaps values of all the doped graphene systems measured from the band structure calculations against different configurations.

It is evident from Fig. 8b that when the pair of Be and S have been substituted for two carbon atoms at different sublattices (non-equivalent) sites with varying relative distances in between the dopants, the electronic band-gaps of the resulting doped systems have been tuned from 0.42 eV through 0.58 eV aside from OB and OC isomer. However, an energy band-gap tuning within an interval of 0.17–0.34 eV has been achieved when the pair of Be and S has been substituted for a pair carbon atoms which are at varying distances apart but at the same sublattices of graphene. This results show that sublattices symmetry is subverted for co-doping at non-equivalent sites but tends to be preserved when is done at equivalent sublattice positions. Therefore Be–S co-doping of graphene creates a finite band-gap in the material. Moreover, the gap can be tuned if the effect of doping sites and dopant–dopant proximity is taken into consideration.

Conclusions

The band-gap of graphene was tuned by co-doping it with Be and S. The sites occupied by the pair of Be and S in graphene were kept at varying distances from each other. The stability, structural and electronic properties of these systems for different lattice sites were determined. The electronic properties of graphene depend on the symmetry, and thus the sites of the impurities play an indispensable role in the modulation of the energy gap. The electronic character changes from semi-metal to semiconductor upon doping and, at the same time the flat structure of the graphene sheet was preserved. The isomers generated by choosing different doping sites differ notably in bond length, band gap and stability. Two set of configurations were considered: one involves Be–S co-doping at equivalent positions while the other configuration deals with substitution of the impurities at non-equivalent sites of graphene sheets. Besides the OB and OC isomer, all the isomers that lent credence to the latter configuration enhanced a band-gap opening but are less stable than the isomers of the Be–S co-doped graphene sheets that shared the platform of the former configuration. As expected, the Fermi level of all the isomers of Be–S co-doped graphene considered in this study lied right in the middle of the band-gap. This shows that Be–S co-doped graphene has no metallic character as the case might be in sulphur (heavily) doped graphene. The result of cohesive and formation energies shows that when Be and S are substituted as nearest neighbour such a conformation gives the most stable configuration than any other one. Due to the differing electronegativity between a carbon atom and the impurities, electrons accumulation occurred on the carbon atoms closest to Be and S in all the doped systems investigated. Our result presents an avenue to tune the electronic band gaps of graphene from 0 eV through 0.58 eV so that the material can be used in electronic applications with operations within this energy range.

Acknowledgements

The author would specifically like to acknowledge Mr Emmanuel Igumbor for his invaluable inputs made to this work. The author would also like to thank Prof. M. E. Casida, and Dr O. Akin-ojo for their suggestions after reading through the manuscript. This work is based on the research supported by the South African Research Chairs Initiative of the Department of Science and Technology and National Research Foundation (NRF) of South Africa (Grant No. 97994). O. Okikiola acknowledges the financial support from NRF and the University of Pretoria for his PhD studies.

References

1. P. Wallace, Phys. Rev., 1947, 71, 622–634.
2. R. Saito, G. Dresselhaus and M. Dresselhaus, Physical properties of carbon nanotubes, 1998.
3. A. K. Geim and K. S. Novoselov, Nat. Mater., 2007, 6, 183–191.
4. A. H. Castro Neto, N. M. R. Peres, K. S. Novoselov, A. K. Geim, A. H. C. Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov and A. K. Geim, Rev. Mod. Phys., 2009, 81, 109–162.
5. N. D. Mermin, Phys. Rev., 1968, 176, 250–254.
6. 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.
7. K. S. A. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. K. I. V. Grigorieva, S. V. Dubonos, A. A. Firsov, M. I. Katsnelson and I. V. Grigorieva, Nature, 2005, 438, 197–200.
8. A. K. Geim, Science, 2009, 324, 1530–1534.
9. E. Castro, H. Ochoa, and M. Katsnelson, R. V. Gorbachev, D. C. Elias, K. S. Novoselov, A. K. Geim and F. Guinea, Phys. Rev. Lett., 2010, 105, 1008.2522.
10. P. A. Denis, Chem. Phys. Lett., 2010, 492, 251–257.
11. J. Zhou, M. Wu, X. Zhou and Q. Sun, Appl. Phys. Lett., 2009, 95, 103108.
12. X. Tian, J. Xu and X. Wang, J. Phys. Chem. B, 2010, 114, 11377–11381.
13. R. Balog, B. Jørgensen, L. Nilsson, M. Andersen, E. Rienks, M. Bianchi, M. Fanetti, E. Lægsgaard, A. Baraldi, S. Lizzit, Z. Sljivancanin, F. Besenbacher, B. Hammer, T. G. Pedersen, P. Hofmann and L. Hornekær, Nat. Mater., 2010, 9, 315–319.
14. L. L. Zhang, X. Zhao, H. Ji, M. D. Stoller, L. Lai, S. Murali, S. Mcdonnell, B. Cleveger, R. M. Wallace and R. S. Ruoff, Energy Environ. Sci., 2012, 5, 9618.
15. P. A. Denis, R. Faccio and A. W. Mombru, ChemPhysChem, 2009, 10, 715–722.
16. J. Dai, J. Yuan and P. Giannozzi, Appl. Phys. Lett., 2009, 95, 232105.
17. I. Zanella, S. Guerini, S. B. Fagan, J. Mendes Filho and A. G. Souza Filho, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 77, 073404.
18. P. Shinde and V. Kumar, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 125401.
19. X. Deng, Y. Wu, J. Dai, D. Kang and D. Zhang, Phys. Lett. A, 2011, 375, 3890–3894.
20. M. Escaño, T. Nguyen and H. Kasai, Chem. Phys. Lett., 2011, 515, 85–90.
21. H. Wang, Q. Wang, Y. Cheng, K. Li, Y. B. Yao, Q. Zhang, C. Z. Dong, P. Wang, U. Schwingenschlogl, W. Yang and X. X. Zhang, Nano Lett., 2012, 12, 141–144.
22. E. Cruz-Silva, Z. Barnett, B. Sumpter and V. Meunier, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 1–9.
23. A. A. Avetisyan, B. Partoens and F. M. Peeters, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 035421.
24. K. Mak, C. Lui, J. Shan and T. Heinz, Phys. Rev. Lett., 2009, 102, 256405.
25. P. Shemella and S. Nayak, Appl. Phys. Lett., 2009, 94, 2101.
26. X. Y. Peng and R. Ahuja, Nano Lett., 2008, 8, 4464–4468.
27. Y. H. Lu, P. M. He and Y. P. Feng, J. Phys. Chem. C, 2008, 112, 12683–12686.
28. L. A. Ponomarenko, F. Schedin, M. I. Katsnelson, R. Yang, E. W. Hill, K. S. Novoselov and A. K. Geim, Science, 2008, 320, 356–358.
29. B. Trauzettel, D. V. Bulaev, D. Loss and G. Burkard, Nat. Phys., 2007, 3, 192–196.
30. M. Han, B. Özyilmaz, Y. Zhang and P. Kim, Phys. Rev. Lett., 2007, 98, 206805.
31. X. Li, X. Wang, L. Zhang, S. Lee and H. Dai, Science, 2008, 319, 1229–1232.
32. J. Bai, X. Zhong, S. Jiang, Y. Huang and X. Duan, Nat. Nanotechnol., 2010, 5, 190–194.
33. S. Ullah, A. Hussain, W. Syed, M. A. Saqlain, I. Ahmad, O. Leenaerts and A. Karim, RSC Adv., 2015, 5, 55762–55773.
34. J. Cai, P. Ruffieux, R. Jaafar, M. Bieri, T. Braun, S. Blankenburg, M. Muoth, A. P. Seitsonen, M. Saleh, X. Feng, K. Müllen and R. Fasel, Nature, 2010, 466, 470–473.
35. S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First, W. A. de Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto and A. Lanzara, Nat. Mater., 2007, 6, 770–775.
36. B. Huang, Q. Xu and S.-H. Wei, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 155406.
37. R. Quhe, J. Zheng, G. Luo, Q. Liu, R. Qin, J. Zhou, D. Yu, S. Nagase, W.-N. Mei, Z. Gao and J. Lu, NPG Asia Mater., 2012, 4, e6.
38. P. Rani and V. K. Jindal, RSC Adv., 2013, 3, 802–812.
39. L. Ferrighi, M. Trioni and C. Di Valentin, J. Phys. Chem. C, 2015, 119, 6056–6064.
40. H. Lee, B. Huang, W. Duan and J. Ihm, J. Appl. Phys., 2010, 107, 084304.
41. Y. Ferro, A. Allouche and C. Linsmeier, J. Appl. Phys., 2013, 113, 213514.
42. P. A. Denis, Comput. Mater. Sci., 2013, 67, 203–206.
43. C. Huang, L. Han, L. Wu, R. Su, J. Chen and P. Lu, Eur. Phys. J. B, 2015, 88, 147.
44. G. K. Gueorguiev, Z. Czigány, A. Furlan, S. Stafström and L. Hultman, Chem. Phys. Lett., 2011, 501, 400–403.
45. W. Tang, E. Sanville and G. Henkelman, J. Phys. Condens. Matter, 2009, 21, 084204.
46. E. Sanville, S. D. Kenny, R. Smith and G. Henkelman, J. Comput. Chem., 2007, 28, 899–908.
47. G. Henkelman, A. Arnaldsson and H. Jónsson, Comput. Mater. Sci., 2006, 36, 354–360.
48. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
49. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
50. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15.
51. G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
52. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
53. P. A. Denis and C. P. Huelmo, J. Phys. Chem. C, 2016, 120, 7103–7112.
54. A. Hussain, S. Ullah and M. Farhan, RSC Adv., 2016, 6, 55990–56003.
55. D. R. Cooper, B. D'Anjou, N. Ghattamaneni, B. Harack, M. Hilke, A. Horth, N. Majlis, M. Massicotte, L. Vandsburger, E. Whiteway and V. Yu, ISRN Condens. Matter Phys., 2012, 2012, 1–56.
56. S. M. Kozlov, F. Viñes and A. Görling, Adv. Mater., 2011, 23, 2638–2643.
57. F. Schwierz, Nat. Nanotechnol., 2010, 5, 487–496.
58. M. Hasegawa and K. Nishidate, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 205431.
59. P. Nath, S. Chowdhury, D. Sanyal and D. Jana, Carbon, 2014, 73, 275–282.
60. P. A. Denis, C. Huelmo and F. Iribarne, Comput. Theor. Chem., 2014, 1049, 13–19.
61. H. L. Poh, P. Šimek, Z. Sofer and M. Pumera, ACS Nano, 2013, 7, 5262–5272.
62. S. Yang, L. Zhi, K. Tang, X. Feng, J. Maier and K. Müllen, Adv. Funct. Mater., 2012, 22, 3634–3640.

---

Footnote

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

---