Structural evolution of BCN systems from graphene oxide towards electrocatalytically active atomic layers

Sumit Bawaria, Kshama Sharmaa, Golap Kalitab, P. K. Madhua, Tharangattu N. Narayanan*a and Jagannath Mondal*a
aTata Institute of Fundamental Research – Hyderabad, Sy. No. 36/P Serilingampally Mandal, Gopanapally Village, Hyderabad-500046, India. E-mail:;;
bDepartment of Physical Science and Engineering, Nagoya Institute of Technology, Gokiso-cho, Showa-ku, Nagoya 466-8555, Japan

Received 7th April 2020 , Accepted 30th May 2020

First published on 2nd June 2020

The boron–carbon–nitrogen (B–C–N) ternary phase diagram is rich with several stable phases. But experimentally exploring the vast array of plausible phases in atomically thin layers is quite challenging. BCN structures vary in nature with the method of preparation, which contributes to differences in their physical properties, particularly while studying photo/electrocatalytic applications. Here, we present a molecular dynamics (MD) approach to study the structural evolution of BCN systems, using the reactive force field (ReaxFF). ReaxFF-based MD simulations allow observation of bond breaking/making in BCN systems, in real time. This helps unravel physically stable structures for B,N-singly and co-doped GO. These predicted structures are subsequently realized in experiments, and characterised via various spectroscopic techniques. We further study the electrochemical activity of these structures for the hydrogen evolution reaction. Specifically, in B,N co-doped systems, we find that the electrocatalytic activity changes considerably if we consider separated doping or domain doping. Together, by combining simulation and experiments, this work shows a structure–property relationship of BCN structures, which ultimately governs their electrocatalytic activity towards the hydrogen evolution reaction.


Heteroatom (B, N, S, P, F, etc.) doped graphene has received tremendous attention in the recent past, primarily due to its electronic and catalytic applications.1,2 Amongst these functional graphene derivatives, B and N individually3,4 and co-doped systems2,5–11 have been highly investigated due to their ‘p’ and ‘n’ dopant nature, respectively. Co-doped systems have received special attention due to the synergistic effect of the heteroatoms and their domains in the parent carbon matrix leading to interesting optical8 as well as spin polarization effects.12 These effects were utilised for photocatalytic applications of BCN structures as well as selective electrochemical reduction of oxygen to peroxide or water.6–8 A few theoretical studies were conducted in the recent past to unravel the role of dopants, co-dopants, and their synergism,13 and they conclude that both electronic and geometric factors can affect the electrocatalytic activity of the doped graphene structures. Although different BCN structures have been experimentally explored towards catalysis, contrasting reports exist. For example, while some studies show high catalytic activity,9,10 others have reported B,N co-doped graphene as a poor catalyst.11 These conflicting reports could be due to microscopic structural variations in the resultant doped structures which arise due to differences in the synthesis protocols and the nature of the precursors.3–12 Though the structural evolution of reduced graphene oxide (rGO) from graphene oxide (GO) was studied previously,14 the structural evolution of atomically thin BCN structures still remains to be explored, thereby lacking mechanistic insight towards their structure–electrocatalytic property correlations.

The electrocatalytic hydrogen evolution reaction (HER) is identified as the cleanest approach towards the production of hydrogen (H2(g)), which is important for ammonia synthesis and also has potential as a high energy density fuel.15 Though the BCN ternary phase diagram is rich with many stable phases of BCN structures,16 the experimental synthesis of these BCN structures via graphene doping is found to form segregated domains of BN and carbon (graphene).17 Out of the different phases of BN, hexagonal boron nitride (hBN) is the only one with a layered nature, and hence atomic layers of BCN mainly contain graphene domains with in-plane connected patchy hBN islands.17 Further, the amount of oxygen can also affect the structural evolution of such BCN structures during thermal annealing (reduction), as has been observed for GO.18 These structural variations in BCN may affect their resultant catalytic performance, from the HER to the oxygen reduction reaction (ORR), where in the ORR both 2e (peroxide) and direct 4e processes are reported for BCN structures.6,7,19

Some of the current authors have previously studied the origin of electrocatalytic and photo-electrocatalytic activity of graphene-hBN van der Waals heterostructures.20–22 In these studies, the roles of the van der Waals shadow effect and possible charge transfer among the individual layers are unveiled as the origin of the HER catalytic activity of such vdW structures.20 Though many studies focus on the ORR activity of graphene-hBN structures, only a few focus on the HER activity of such in-plane structures.11 For example, Jiao et al. have demonstrated the origin of the HER activity in heteroatom doped graphene (B, N, S, and P doped structures) using density functional theory (DFT) based calculations.23 Their studies propose an experimentally achievable two step method towards graphene based HER catalysts, having comparable performance with benchmark catalysts via intrinsic electronic property modulation and an increased number of active sites for an identified doped configuration.23 But studies on such ideal structures using DFT are still far away from the experimental structures, for example, where most of the large scale BCN samples were produced via thermal reduction of GO in the presence of B and N sources.9–11 Under controlled conditions, it is possible to form specific BCN configurations between graphene and hBN and the domain size can also be tuned; however, these processes are not scalable.24,25 The presence of oxygen functionalities can affect the resultant microscopic structure of BCN, and hence studies on the origin of the HER on real structures are highly demanding.

Studying these complex systems using DFT is very computationally expensive, and classical molecular dynamics (MD) simulations cannot predict new structures, as no bond breaking/formation is allowed. To this end, in this work, we use MD with a reactive force field (ReaxFF),26 which provides a DFT-based bond-order parameter that enables us to study bond breaking/formation on long timescales. Using this approach, unbiased structural evolution of different atomically thin BCN structures is carried out. As mentioned before, Bagri et al. have employed a reactive MD approach to ascertain whether MD can reproduce the functional groups seen in experiments and perhaps could point out unknown functionalities.14 Here, we use a similar strategy to uncover the possible functional groups on B,N individually and co-doped systems formed via the annealing of GO in the presence of B and N precursors. Various spectroscopic (micro-Raman, FTIR, NMR, and XPS) and microscopic (TEM) techniques are used to correlate the theoretical findings with the experimentally developed samples. Finally, the effect of structural variations in BCN structures towards their electrocatalytic HER activity is subjected to investigation and the results are correlated with the DFT-based charge-transfer studies.


Reactive MD simulations

A system consisting of a periodic ∼4 × 4 nm graphene sheet was constructed, consisting of 680 atoms. 20% and 33% O/C ratio structures were made by adding alcohol (C–OH) and epoxy (C–O–C) groups in close proximity as suggested by Cai et al.27 B/N doped systems consider ∼5% of the carbons replaced by B/N, and in the case of B,N doped systems 5% B and 5% N were introduced. BN-GO contained a ∼1 × 2 nm domain of BN, consisting of 125 atoms. All MD simulations were performed using the open source computation program LAMMPS.28,29 To account for bonding interactions, the ReaxFF force-field,28 as developed by Paupitz et al.,30 was used for reactive MD simulations, which consisted of bonding parameters between C, B, N, O, and H. MD simulations were run using an NVT ensemble in a box of volume 4.17 × 4.25 × 2.0 nm3, and the Berendsen thermostat was used to modulate the temperature. The annealing process consisted of a 250 fs initialisation to 1500 K, 200–500 ps annealing at 1500 K, and a final 1.25 ps quenching to 300 K.

Experimental preparation of doped GO

GO was prepared using the modified Hummers' method.31 Doping of GO was performed using a chemical vapor deposition (CVD) chamber assisted high temperature annealing method, where the sample was kept in a ceramic boat, as shown in Fig. S1 (ESI). The sample was annealed at 800 °C for 30 minutes in an N2(g) atmosphere. For B-GO and N-GO, GO was mixed with boric acid and melamine in a 1[thin space (1/6-em)]:[thin space (1/6-em)]5 w/w ratio, respectively. For B,N-GO, as prepared B-GO is mixed with melamine in a 1[thin space (1/6-em)]:[thin space (1/6-em)]5 w/w ratio, followed by high temperature annealing at 800 °C. For the co-doped BN-GO development, boric acid, melamine, and GO were kept sequentially such that evaporating boric acid and melamine react with each other and hence the B and N sources will reach the GO surface simultaneously.


Fourier transform infrared (FTIR) spectra were recorded in transmittance mode using a Bruker (model: α) spectrometer. A Renishaw inVia micro-Raman spectrometer with a 532 nm laser was used (10 s acquisition at 1% intensity) for micro-Raman analyses. For XPS, an ESCA Plus spectrometer (Omicron Nanotech Ltd) was used, with a monochromated Mg Kα source and 1.2 keV excitation. For TEM, the sample structure analysis and elemental mapping were carried out using a FE-TEM JEOL:JEM-2100F, operated at 200 kV.


All solid-state NMR spectra were acquired on a Bruker-Avance 500 (11.7 T) NMR spectrometer using a 4 mm MAS probe with samples fully packed inside zirconium oxide rotors. The MAS spinning frequency was 8.333 kHz for 13C and 11B single-pulse (Bloch decay) experiments. For the single-pulse 13C NMR experiment, 6900 transients were averaged with a recycle delay of 10 s with an acquisition time of 5.2 ms. All 13C spectra were referenced externally to adamantane for 13C (δref = 38.4 ppm for –CH2 of adamantane). For all 11B single-pulse experiments, the radio frequency (RF) amplitude was such that the pulse length was 2.0 μs, and a recycle delay time of 0.5 s was employed with an acquisition time of 5 ms. The chemical shifts for 11B were referenced using 15% BF3.OEt2 in CDCl3 (δref = 0 ppm). Data were processed using Topspin 3.5 with an exponential line broadening of 30 Hz.

Electrochemical HER studies

All the electrochemical studies were performed using a Biologic SP-300 potentiostat. A three-electrode system consisting of a counter electrode (platinum/graphite rod), reference electrode (Ag/AgCl), and working electrode (glassy carbon electrode (GCE)) was used, in 0.5 M H2SO4. The GCE was modified with different materials via drop-coating. An activation cyclic voltammogram (CV) of 20 cycles was run at 50 mV s−1, preceding an LSV at 20 mV s−1.


Structures that are analogues of the MD annealed systems were studied using the SIESTA 4.0 open source code.32 Pseudopotentials were generated using the norm-conserving Troullier–Martins scheme, with the GGA-PBE exchange correlation functional.33,34 The energy cut-off for the real space grid was set at 500 Ry. Variable-cell relaxation was carried out using the Broyden method, till the maximum force on each atom was less than 0.00001 Hartree Bohr−1.35 A 4 × 4 × 1 k-point sampling of the Brillouin zone was carried out using the Monkhorst–Pack scheme.36

Results and discussion

ReaxFF simulations of doped GO and experimental annealing

We begin with the thermal annealing of graphene oxide (GO).14 Two initial O/C (oxygen/carbon) ratios are considered here, namely 20% and 33% O/C, which represent low- and high-oxygen containing GO samples. The initial oxygen functionalities consist entirely of epoxy (C–O–C) and alcohol (C–OH) groups in proximity to each other, as suggested by Cai et al.25 The initial and annealed structures for GO 20% and 33% O/C are shown in Fig. S2 (ESI). Starting from simply epoxy and alcohol groups, the thermal annealing process was able to generate diverse functional groups, like carbonyls, in-plane oxygens (pyran, furan, etc.), lactols, and phenol, by using reactive molecular dynamics simulation (Fig. S3, ESI).14 To simulate the thermal annealing process of GO with a B/N precursor, we considered four doped systems, namely B-doped (B-GO), N-doped (N-GO), B,N separately doped (B,N-GO), and B,N domain doped GO (BN-GO). All systems were annealed at 1500 K, and the final structures are shown in Fig. 1. As the structure is annealed, much of the oxygen species dissociate from the graphitic lattice, in the form of H2O, O2, CO2, etc. Finally, we see different kinds of functionalities with an overall graphitic lattice structure for 20% (Fig. 1(1–4)), and a lot more disorder for 33% (Fig. 1(5–8)) structures.
image file: d0qm00220h-f1.tif
Fig. 1 Final structures obtained in reactive MD simulations, after annealing at 1500 K, for B-GO (1,5), N-GO (2,6), B,N-GO (3,7) and BN-GO (4,8) at 20% (1–4) and 33% (5–8) O/C ratios.

For experimental realisation, we synthesised BCN structures by annealing in a CVD (chemical vapour deposition) chamber. To experimentally dope B and/or N in GO, GO is annealed with boric acid/melamine at 800 °C under an N2 atmosphere (Fig. S1, ESI), as explained in the methods. Thermal reduction of GO is apparent by the decrease in the OH (∼3430 cm−1) and C[double bond, length as m-dash]O (∼1725 cm−1) peaks in the FTIR spectra (Fig. S4, ESI).3 For B-GO and N-GO, a right shift in the C–C bond stretching is clearly seen (∼1630 to ∼1565 cm−1), which indicates doping of B/N.3 FTIR spectra of BN co-doped samples, synthesised either sequentially (B,N-GO) or simultaneously (BN-GO), are also shown in Fig. S2 (ESI). For B,N-GO and BN-GO, a clear BN stretching peak (∼1400 cm−1) is observed. This indicates the possibility of in-plane B and N doping in rGO lattices. However, there can be multiple possible functionalities in each of these samples. The formation of atomically thin structures of BCN is evident from the TEM analyses (Fig. S5, ESI).

Table 1 compares the B, C, N, and O percentages obtained using XPS (Fig. S6, ESI) for the annealed samples and of the simulated structures. We find reasonably good correlation between experiments and simulations, with the doping percentage in B-GO and N-GO being close to 2–3%, compared to the ∼4% considered in the simulations. The oxygen percentage also matches remarkably well for the 20% O/C simulations. B,N-GO shows a much higher oxygen concentration due to the presence of boron oxide (B2O3) and nitrogen oxide (NOx) side products, as seen in high-resolution XPS (Fig. 3A and B). BN-GO shows the highest B and N ratios at ∼5.5%, which may be due to the higher stability of BN domains. The overall elemental percentage and changes in oxygen percentage found in XPS are well reflected in the 20% MD structures. To compare the different functionalities in simulated and experimentally annealed samples, further analyses are performed.

Table 1 Comparison between B, C, N, and O percentages between experimental (XPS) and MD simulated, annealed samples
System % XPS % MD (20%)
B 2.46 4.10
C 83.23 82.95
O 14.31 12.95
N 2.84 4.14
C 84.50 86.11
O 12.66 9.75
B 4.85 3.75
N 3.92 4.14
C 66.83 79.56
O 24.40 12.55
B 5.48 8.64
N 5.58 8.77
C 78.48 76.74
O 10.46 5.85

Defect density and oxygen functionalities

The percentage of defects in all simulated structures is calculated by counting the nearest neighbours within 0.18 nm of each carbon atom. A carbon with three carbon neighbours is considered as an ideal sp2 type carbon. Fig. 2A shows the percentage of non-sp2 type carbons, which can be considered as the defect percentage. As expected, the 33% O/C structures are much more defective than the 20% structures. Of all four 20% O/C structures, BN-GO and N-GO are the ones with the least defects, followed by BGO and B,N-GO as the most defective structures. This correlates well with micro-Raman spectra (Fig. 2C), which show a lower ID/IG (intensity of D peak/intensity of G peak) ratio (0.6) for NGO and BN-GO, compared to the ID/IG ratio (0.8) for B-GO and B,N-GO. In addition to this, the Raman spectrum of BN-GO shows a peak at ∼1370 cm−1 (Fig. 2C) that can be ascribed to BN domain vibrations, which is shifted from 1364 cm−1 for pristine hBN. No B–N vibration peak (∼1364 cm−1) is seen in B,N-GO. The MD-derived defect percentage for most 33% samples closely follows the 20% trend, except N-GO 33%, which is much more defective than expected. This difference can be ascribed to the fact that nitrogen does not bond easily with oxygen functionalities, in contrast to boron, which in turn leads to repulsion between N and O, and more structural defects.
image file: d0qm00220h-f2.tif
Fig. 2 (A) % Non-graphitic (non-sp2) carbons or the defect percentage for B-GO, N-GO, B,N-GO and BN-GO in 20% and 33% O/C ratios. (B) Amount of different oxygen functionalities for all the samples. (C) Raman spectra for all the samples showing the D peak (∼1350 cm−1), G peak (∼1605 cm−1), and the BN peak (∼1370 cm−1), which are highlighted. (D) Deconvoluted C 1s XPS spectra for all the four samples.

The number of different functional groups (C[double bond, length as m-dash]O, C–OH, C–O–C, and B–O–C) that oxygen forms is shown in Fig. 2B. In the basal plane of all samples, epoxy groups and carboxyl groups dominate, with intermittent in-plane bonded oxygens. Going from 20% to 33%, a major rise in OH groups is seen due to the ‘ripping’ of sheets. These oxygen functionalities can be compared to the solid-state NMR spectrum of GO (Fig. S7, ESI), which shows a high presence of C–O–C (61 ppm) (in-plane and epoxy) and C–OH (70 ppm) (alcohol) functional groups, amidst C–C peaks for both reduced (112 ppm) and oxidised (130 ppm) regions.17 C[double bond, length as m-dash]O (194 ppm) (carbonyl), O–C[double bond, length as m-dash]O (169 ppm) (carboxyl), and O–C–O (101 ppm) (lactol) functional groups are also present in lower concentrations. All these functionalities are well represented in the simulated structure of GO (Fig. S3, ESI), and also all doped structures. However, since we only consider the basal plane in our simulations, –OH groups are under-represented, which are found in abundance at the edges. Also, in the case of BGO and B,N-GO, most of the boron species are directly bonded to oxygen (B–O–C). High resolution XPS spectra for carbon show the presence of C[double bond, length as m-dash]C, C–O, C[double bond, length as m-dash]O, and O–C[double bond, length as m-dash]O species in all samples (Fig. 2D), which compares with the C[double bond, length as m-dash]O, C–OH, and C–O–C functionalities found in MD. Apart from oxygen functional groups, the presence of C–B in B-GO and BN-GO, and C–N in N-GO, B,N-GO, and BN-GO is also confirmed through XPS. The diverse functional groups formed by boron and nitrogen are studied in more detail.

The nature of boron and nitrogen functionalities

For all simulated structures, the initial structure has only B/N point defects, with oxygen species in proximity. As the structure was annealed, different functional groups evolved. The XPS spectra for B-GO (Fig. 3A) can be deconvoluted into three individual peaks of B–3C (191.27 eV), B–(2C)–O (192.86 eV), and B–(C)–2O (194.67 eV), which can be directly compared to the simulated structure. The percentage of B–3C type (∼12.64%) in XPS analysis is in excellent agreement with that obtained from ReaxFF MD simulations (∼10.31%). Some representative snapshots of boron functionalities are shown in Fig. 3A(i–iv). Since boron has a high tendency to bind with oxygen, B–(O)–C type functionalities dominate. In the case of N-GO, the XPS spectrum can be deconvoluted into graphitic N–3C (401.51 eV), pyridinic N–2C (398.16 eV), and pyrrolic N–2C (399.74 eV) species. Representative snap shots (Fig. 3B) show a graphitic (i), pyridinic (ii), cyano (iii), and pyrrolic (iv) nitrogen in N-GO. Our MD simulation, which starts from completely graphitic type nitrogen, yields very few pyrrolic functionalities. And if we compare the N–3C type functionality, the XPS derived (∼26.08%) and MD simulated (∼32.45%) results are not very far from each other. As shown in Fig. 2B, nitrogen does not easily bind to oxygen, and hence very few N–O type functionalities are seen in the simulation. In the case of B,N-GO, we consider a system that is simply a mix of B-GO and N-GO type functionalities. However, the B 1s and N 1s XPS spectra of B,N-GO show a lot more oxygenated functionalities, which explain the high oxygen content (Table 1). In the case of BN-GO, the B 1s and N 1s XPS spectra show the presence of similar functionalities, as in B-GO and N-GO. Moreover, new regions that correspond to B–N bonding are seen to be enhanced, showing BN domain formation.
image file: d0qm00220h-f3.tif
Fig. 3 (A) Deconvoluted B 1s XPS spectra of B-GO, B,N-GO and BN-GO, along with representative boron functional groups (i–iv). (B) Deconvoluted N 1s XPS spectra of NGO, B,N-GO, and BN-GO, along with representative nitrogen functional groups (i–iv). (C) Solid-state 11B NMR spectra for B-GO, B,N-GO, and BN-GO, with B–C (B4C), B–O (H3BO3), and B–N (hBN) peaks highlighted. (D) 13C NMR spectra for all the four samples.

To confirm the local environment around boron, we performed single-pulse 11B solid-state NMR for B-GO, B,N-GO, and BN-GO, along with B4C (B–C), boric acid (B–O), and hBN (B–N) powders as reference samples (Fig. 3C). In B-GO, only B–C (∼0.17 ppm) and B–O (∼10.83 ppm) type functionalities are seen, with very broad peaks. B,N-GO shows a highly resolved (may be ordered) spectrum, and two distinct regions are seen for B–C (∼0.17 ppm) and B–O (∼10.83 ppm), with a slight shoulder for B–N (∼18.11 ppm). Meanwhile for BN-GO, the spectrum is similar to B-GO, with a second peak region due to B–N interactions (∼18.11 ppm). The NMR spectrum suggests considerable B–N formation (BN domains) in BN-GO, in tune with the observations of FTIR, Raman, and XPS analyses. However the NMR spectra of B,N-GO do not show any substantial presence of domains, in tune with fact that no B–N peaks are found in XPS analyses, and hence these studies suggest that the structure has majorly separated B and N dopants, along with B/N-oxide side-products. The 13C NMR spectra for all four samples (Fig. 3D) show only a single peak at 112 ppm, which corresponds to C–C bonding in reduced regions. The only difference is in the line widths (full width at half maximum), with N-GO (368.04 Hz) and BN-GO (351.72 Hz) having narrow line widths, and B-GO (540 Hz) and B,N-GO (499.27 Hz) having a much wider line width. The broadening of the spectra of B-GO and B,N-GO (and also the slight shift in the peak of B,N-GO in comparison to that of BN-GO) is an indication of their disordered nature, which is also evident from Raman and MD studies. Structural variations in all four samples are evident from these studies, which can affect their properties.

Electrocatalytic HER

Next we highlight the impact of these structural variations in the electrochemical HER performance of BCN systems. The HER studies were conducted in acidic medium (pH = 0) and the linear sweep voltammetry (LSV) HER polarisation curves are shown in Fig. 4A. It is apparent from the HER LSV curves that BN-GO is the ‘best catalyst’ out of all the samples in terms of lower overpotential and higher current density. BN-GO has an overpotential of ∼0.2 V vs. RHE. N-GO and B,N-GO perform with a similar overpotential ∼0.4 V, followed by B-GO (0.5 V). The HER overpotentials at a current of 2 mA cm−2 (set current density for comparison of overpotentials) measured here are in close agreement with those measured by others (mainly for N-GO and B,N-GO).11,23 Here, we differentiate codoped systems as B,N-GO11,23 and BN-GO,9,10 depending on whether the B/N sources are used sequentially or simultaneously, respectively. The overpotential for B-GO shows different values in three different studies (including our own).4,23 It is apparent from all these studies that BN-GO is the ‘best’ performing HER catalyst out of these different B/N GO systems, which is likely due to the formation of less defective BN domain containing rGO structures. We find an anti-correlation between the defect density and experimental HER, with the lower ID/IG ratio sample, i.e. BN-GO and N-GO, performing better than B-GO and B,N-GO. This is in contrast to an earlier report,37 suggesting that not all types of defects are advantageous towards catalytic processes.
image file: d0qm00220h-f4.tif
Fig. 4 (A) Electrochemical HER LSV polarisation curves for GO, B-GO, N-GO, B,N-GO and BN-GO, performed in 0.5 M H2SO4. (B) Comparison of overpotentials with existing literature.4,9–11,23 (C) ΔGH values for different systems, in B-GO, N-GO and BN-GO. (D) Weighted ΔGH (weighted using MD concentrations), plotted against the experimental overpotential.

DFT calculations on MD analogs

As seen from our reactive MD simulated structures (Fig. 1), the structure of doped GO is much more complicated than single point substitutions of B/N in the graphitic lattice. We see oxygenated boron functionalities in B-GO, and pyridinic/pyrrolic nitrogens in N-GO, apart from graphitic substitutions. To account for this diversity in functional groups, we consider all major functionalities found in MD, and perform ΔGH calculations on multiple sites (Fig. S8, ESI). The lowest ΔGH values found for all systems are shown in Fig. 4C, which show a range of values. B-GO consists of B–3C (0.10 eV), B–O–2C (0.81 eV), and B–2O–C (0.81 eV) type sites. In the case of N-GO, we study graphitic (0.31 eV), pyridinic (0.59 eV), and pyrrolic (0.33 eV) sites. When BN-GO is studied using an ideal BN domain inside graphene, the ΔGH values are quite high (>0.9 eV). However, we find from the MD annealed structures for BN-GO (Fig. 1) that the interfaces between the domains of carbon and BN are rather defective. If we consider defective boundaries, we obtain lower ΔGH values for a boron deficiency (B-def) (0.11 eV) and a nitrogen deficiency (N-def) (0.47 eV).

The ΔGH value obtained for each system was weighted according to the concentration of that species found in MD simulations (Table S2, ESI). Fig. 4D shows the weighted ΔGH for different systems correlated with the overpotential of the experimental sample. BN-GO shows the expected lowest weighted ΔGH and the lowest overpotential, followed by N-GO, B,N-GO and B-GO. This difference shows a clear distinction between BN-GO and B,N-GO, which are both B,N-codoped systems. It is apparent from Fig. 4D that the weighted ΔGH value correlates well with the experimental overpotential.

Further, we found that simulating the doping process with boron and nitrogen sources such as borane and melamine (instead of direct boron and nitrogen substitution), respectively, also gives similar doped structures (Fig. S9, ESI), though the amount of dopant is low. For example, in the case of boron, B–O–C type bonding is found to be prevalent, which is very similar to that in BGO obtained via elemental substitution. Though nitrogen initially forms C–NH2 type species when annealing GO with melamine, eventually it is found to be forming graphitic/pyridinic species. Hence, both of these results are in agreement with elemental doping, further validating the ReaxFF based structural variations in the BCN system, as discussed above.


In conclusion, we have simulated annealed structures of B/N doped graphene oxide and substantiated them with FTIR, Raman, XPS, and NMR measurements. All the four structures (B-GO, N-GO, B,N-GO, and BN-GO) are dominated by different functionalities. The structures are further used to understand the different active sites in each sample, and model them using DFT to determine the HER activity (ΔGH). A previous report comparing ΔGH for B-G and N-G shows lower ΔGH for B-G, while N-G is clearly the better catalyst in actual experiments.23 We expect that this discrepancy arises due to the fact that the site chosen for ΔGH calculation is under-represented in the actual sample. In the case of B-GO, although a B–3C type site has a low ΔGH (0.10 eV), B–O–2C and B–2O–C type sites (0.81 eV) dominate, and bring down the catalytic activity. This suggests that a more reduced B-GO type sample, with higher B–3C, will be a better catalyst. An overall good correlation between the DFT-derived weighted ΔGH and experimental overpotential is found for all samples. A clear distinction is found between B,N-GO (separately doped) and BN-GO (domain doped) codoped systems, which suggests a role played by the domain boundary. This study highlights the significance of understanding the atomic structure of a catalyst, which can help explain subtle differences in electrocatalytic activity.

Conflicts of interest

There are no conflicts to declare.


We are thankful for the support of this project from the intramural grants at TIFR Hyderabad from the Department of Atomic Energy (DAE). The authors thank Prof. M. M. Shaijumon and his group at Indian Institute of Science, Education, and Research – Thiruvananthapuram, India for helping to conduct the XPS measurements.


  1. R. Lv and M. Terrones, Towards new graphene materials: doped graphene sheets and nanoribbons, Mater. Lett., 2012, 78, 209–218 CrossRef CAS.
  2. J. Duan, S. Chen, M. Jaroniec and S. Z. Qiao, Heteroatom-Doped Graphene-Based Materials for Energy-Relevant Electrocatalytic Processes, ACS Catal., 2015, 5, 5207–5234 CrossRef CAS.
  3. X. Xu, T. Yuan, Y. Zhou, Y. Li, J. Lu, X. Tian, D. Wang and J. Wang, Facile synthesis of boron and nitrogen-doped graphene as efficient electrocatalyst for the oxygen reduction reaction in alkaline media, Int. J. Hydrogen Energy, 2014, 39, 16043–16052 CrossRef CAS.
  4. B. R. Sathe, X. Zou and T. Asefa, Metal-free B-doped graphene with efficient electrocatalytic activity for hydrogen evolution reaction, Catal. Sci. Technol., 2014, 4, 2023–2030 RSC.
  5. S. Wang, L. Zhang, Z. Xia, A. Roy, D. W. Chang, J. B. Baek and L. Dai, BCN graphene as efficient metal-free electrocatalyst for the oxygen reduction reaction, Angew. Chem., Int. Ed., 2012, 51, 4209–4212 CrossRef CAS PubMed.
  6. Y. Zheng, Y. Jiao, L. Ge, M. Jaroniec and S. Z. Qiao, Two-step boron and nitrogen doping in graphene for enhanced synergistic catalysis, Angew. Chem., Int. Ed., 2013, 52, 3110–3116 CrossRef CAS PubMed.
  7. S. Chen, Z. Chen, S. Siahrostami, D. Higgins, D. Nordlund, D. Sokaras, T. R. Kim, Y. Liu, X. Yan, E. Nilsson, R. Sinclair, J. K. Nørskov, T. F. Jaramillo and Z. Bao, Designing Boron Nitride Islands in Carbon Materials for Efficient Electrochemical Synthesis of Hydrogen Peroxide, J. Am. Chem. Soc., 2018, 140, 7851–7859 CrossRef CAS PubMed.
  8. M. Li, Y. Wang, P. Tang, N. Xie, Y. Zhao, X. Liu, G. Hu, J. Xie, Y. Zhao, J. Tang, T. Zhang and D. Ma, Graphene with Atomic-Level In-Plane Decoration of h-BN Domains for Efficient Photocatalysis, Chem. Mater., 2017, 29, 2769–2776 CrossRef CAS.
  9. M. Chhetri, S. Maitra, H. Chakraborty, U. V. Waghmare and C. N. R. Rao, Superior performance of borocarbonitrides, BxCyNz, as stable, low-cost metal-free electrocatalysts for the hydrogen evolution reaction, Energy Environ. Sci., 2016, 9, 95–101 RSC.
  10. S. Ozden, S. Bawari, S. Vinod, U. Martinez, S. Susarla, C. Narvaez, J. Joyner, C. S. Tiwary, T. N. Narayanan and P. M. Ajayan, Interface and defect engineering of hybrid nanostructures toward an efficient HER catalyst, Nanoscale, 2019, 11, 12489–12496 RSC.
  11. Y. Zheng, Y. Jiao, L. H. Li, T. Xing, Y. Chen, M. Jaroniec and S. Z. Qiao, Toward design of synergistically active carbon-based catalysts for electrocatalytic hydrogen evolution, ACS Nano, 2014, 8, 5290–5296 CrossRef CAS PubMed.
  12. M. Fan, J. Wu, J. Yuan, L. Deng, N. Zhong, L. He, J. Cui, Z. Wang, S. K. Behera, C. Zhang, J. Lai, B. M. I. Jawdat, R. Vajtai, P. Deb, Y. Huang, J. Qian, J. Yang, J. M. Tour, J. Lou, C. W. Chu, D. Sun and P. M. Ajayan, Doping Nanoscale Graphene Domains Improves Magnetism in Hexagonal Boron Nitride. Advanced Materials, Adv. Mater., 2019, 31, 1–9 Search PubMed.
  13. J. F. Chen, Y. Mao, H. F. Wang and P. Hu, Theoretical Study of Heteroatom Doping in Tuning the Catalytic Activity of Graphene for Triiodide Reduction, ACS Catal., 2016, 6, 6804–6813 CrossRef CAS.
  14. A. Bagri, C. Mattevi, M. Acik, Y. J. Chabal, M. Chhowalla and V. B. Shenoy, Structural evolution during the reduction of chemically derived graphene oxide, Nat. Chem., 2010, 2, 581–587 CrossRef CAS PubMed.
  15. J. A. Turner, Sustainable Hydrogen Production, Science, 2004, 305, 972–974 CrossRef CAS PubMed.
  16. L. Song, Z. Liu, A. L. M. Reddy, N. T. Narayanan, J. Taha-Tijerina, J. Peng, G. Gao, J. Lou, R. Vajtai and P. M. Ajayan, Binary and ternary atomic layers built from carbon, boron, and nitrogen, Adv. Mater., 2012, 24, 4878–4895 CrossRef CAS PubMed.
  17. L. Ci, L. Song, C. Jin, D. Jariwala, D. Wu, Y. Li, A. Srivastava, Z. F. Wang, K. Storr, L. Balicas, F. Liu and P. M. Ajayan, Atomic layers of hybridized boron nitride and graphene domains, Nat. Mater., 2010, 9, 430–435 CrossRef CAS PubMed.
  18. W. Gao, L. B. Alemany, L. Ci and P. M. Ajayan, New insights into the structure and reduction of graphite oxide, Nat. Chem., 2009, 1, 403–408 CrossRef CAS PubMed.
  19. P. K. Rastogi, K. R. Sahoo, P. Thakur, R. Sharma, S. Bawari, R. Podila and T. N. Narayanan, Graphene-hBN non-van der Waals vertical heterostructures for four-electron oxygen reduction reaction, Phys. Chem. Chem. Phys., 2019, 21, 3942–3953 RSC.
  20. S. Bawari, N. M. Kaley, S. Pal, T. V. Vineesh, S. Ghosh, J. Mondal and T. N. Narayanan, On the hydrogen evolution reaction activity of graphene-hBN van der Waals heterostructures, Phys. Chem. Chem. Phys., 2018, 20, 15007–15014 RSC.
  21. S. Bawari, T. N. Narayanan and J. Mondal, Atomistic Elucidation of Sorption Processes in Hydrogen Evolution Reaction on a van der Waals Heterostructure, J. Phys. Chem. C, 2018, 122, 10034–10041 CrossRef CAS.
  22. S. Bawari, S. Pal, S. Pal, J. Mondal and T. N. Narayanan, Enhanced Photo-Electrocatalytic Hydrogen Generation in Graphene/hBN van der Waals Structures, J. Phys. Chem. C, 2019, 123, 17249–17254 CrossRef CAS.
  23. Y. Jiao, Y. Zheng, K. Davey and S. Z. Qiao, Activity origin and catalyst design principles for electrocatalytic hydrogen evolution on heteroatom-doped graphene, Nat. Energy, 2016, 1, 16130 CrossRef CAS.
  24. Y. Gong, G. Shi, Z. Zhang, W. Zhou, J. Jung, W. Gao, L. Ma, Y. Yang, S. Yang, G. You, R. Vajtai, Q. Xu, A. H. Macdonald, B. I. Yakobson, J. Lou, Z. Liu and P. M. Ajayan, Direct chemical conversion of graphene to boron- and nitrogen- and carbon-containing atomic layers, Nat. Commun., 2014, 5, 3193 CrossRef PubMed.
  25. Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Lei, X. Yang, J. Zhang, J. Yu, K. P. Hackenberg, A. Babakhani, J. C. Idrobo, R. Vajtai, J. Lou and P. M. Ajayan, In-plane heterostructures of graphene and hexagonal boron nitride with controlled domain sizes, Nat. Nanotechnol., 2013, 8, 119–124 CrossRef CAS PubMed.
  26. T. P. Senftle, S. Hong, M. M. Islam, S. B. Kylasa, Y. Zheng, Y. K. Shin, C. Junkermeier, R. Engel-Herbert, M. J. Janik, H. M. Aktulga, T. Verstraelen, A. Grama and A. C. T. van Duin, The ReaxFF reactive force-field: development, applications and future directions, npj Comput. Mater., 2016, 2, 15011 CrossRef CAS.
  27. W. Cai, R. D. Piner, F. J. Stadermann, S. Park, M. A. Shaibat, Y. Ishii, D. Yang, A. Velamakanni, J. A. Sung, M. Stoller, J. An, D. Chen and R. S. Ruoff, Synthesis and solid-state NMR structural characterization of 13C-labeled graphite oxide, Science, 2008, 321, 1815–1817 CrossRef CAS PubMed.
  28. S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, J. Comput. Phys., 1995, 117, 1–19 CrossRef CAS.
  29. H. M. Aktulga, J. C. Fogarty, S. A. Pandit and A. Y. Grama, Parallel reactive molecular dynamics: numerical methods and algorithmic techniques, Parallel Computing, 2012, 38, 245–259 CrossRef.
  30. R. Paupitz, C. E. Junkermeier, A. C. T. van Duin and P. S. Branicio, Fullerenes generated from porous structures, Phys. Chem. Chem. Phys., 2014, 16, 25515–25522 RSC.
  31. D. C. Marcano, D. V. Kosynkin, J. M. Berlin, A. Sinitskii, Z. Z. Sun, A. Slesarev, L. B. Alemany, W. Lu and J. M. Tour, Improved Synthesis of Graphene Oxide, ACS Nano, 2010, 4, 4806–4814 CrossRef CAS PubMed.
  32. J. M. Soler, E. Artacho, J. D. Gale, A. García, J. Junquera, P. Ordejón and D. Sánchez-Portal, The SIESTA method for ab initio order-N materials simulation, J. Phys.: Condens. Matter, 2002, 14, 2745 CrossRef CAS.
  33. N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1991, 43, 1993 CrossRef CAS PubMed.
  34. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865 CrossRef CAS PubMed.
  35. D. D. Johnson, Modified Broyden's method for accelerating convergence in self-consistent calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 38, 12807 CrossRef PubMed.
  36. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188 CrossRef.
  37. J. H. Zhong, J. Zhang, X. Jin, J. Y. Liu, Q. Li, M. H. Li, W. Cai, D. Y. Wu, D. Zhan and B. Ren, Quantitative correlation between defect density and heterogeneous electron transfer rate of single layer graphene, J. Am. Chem. Soc., 2014, 136, 16609–16617 CrossRef CAS PubMed.


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

This journal is © the Partner Organisations 2020