Reassessing structural models of graphitic carbon nitride for reliable photocatalytic predictions

Abstract

Among emerging photocatalysts, graphitic carbon nitride (g-C₃N₄) has attracted considerable attention as a metal-free, visible-light-active semiconductor for sustainable hydrogen production and environmental remediations. Here we demonstrate that many of the obvious models for the layer structure and preferred stacking arrangements of g-C₃N₄ lead to incorrect predictions of photocatalytic properties. We present a systematic structural exploration of monolayer and bulk geometries of heptazine-based g-C₃N₄ using first principles density functional theory (DFT) to assess their stability and elucidate how these configurations influence their electronic properties. Across all systems investigated, we show that using a buckled heptazine structure—rather than the conventional planar geometry—represents the true energetic ground-state structure. A buckled structure is also found by reverse Monte Carlo analysis of experimental data, also reported here. Further analysis of buckling and different stacking registries in g-C₃N₄ layered structures reveals newly identified low-energy corrugated stackings (with P1 symmetry) that give rise to a broad range of electronic band gaps, several of which closely match experimentally reported values of ∼2.7 eV. Moreover, introducing non-metal P and metal Ni dopants at identical lattice sites across 2D-planar/buckled and 3D-corrugated hosts shows distinct electronic responses. Our findings demonstrate that accurate structural modelling, including buckling and stacking registry, is crucial for reliably capturing the stability and tunable electronic properties of g-C₃N₄.

---

1 Introduction

The urgent demand for sustainable energy generation and environmental protection has accelerated interest in solar-driven photocatalysis as a practical route to mitigate these global challenges.¹ Semiconductor-based photocatalysts, capable of harvesting sunlight to drive chemical reactions such as water splitting for hydrogen production, carbon dioxide reduction, and degradation of organic pollutants, provide a clean and scalable approach towards a clean energy future.^2,3 Among the wide range of semiconductor photocatalysts investigated, graphitic carbon nitride (g-C₃N₄) has emerged as a promising two-dimensional (2D) metal-free polymeric photocatalyst due to its suitable bandgap for visible light activity, chemical stability, and ease of synthesis, thereby making it highly relevant for various sustainable applications.^4,5

Graphitic carbon nitride exists in two structural frameworks, comprising either s-triazine (C₃N₃) units or tri-s-triazine (C₆N₇, also known as heptazine) units. Both structures consist of fused 6-membered rings with sp² hybridized C and N atoms, forming delocalized π-conjugated polymeric networks. Moreover, as revealed by X-ray diffraction (XRD), g-C₃N₄ is a hexagonal layered structure with two characteristic peaks at 2θ of approximately 13.04° and 27.6°, attributed to the 〈101̄0〉 in-plane periodicity and 〈0002〉 interlayer stacking governed by weak van der Waals interactions, respectively. While s-triazine and heptazine share a layered graphitic arrangement, heptazine units, consisting of three fused triazine rings, are energetically more favourable.^6–8 Zhang et al.⁹ report the structural evolution from triazine to heptazine during thermal polymerization and highlights new theoretical insights for synthesizing heptazine-based g-C₃N₄ at lower temperatures.

Unlike graphene, which consists mainly of carbon atoms, polymeric g-C₃N₄ contains both carbon and nitrogen atoms, resulting in increased structural complexity. Consequently, uncertainty regarding the precise structure of polymeric g-C₃N₄ and the identification of its true ground-state structure continues to present significant challenges for both experimental and theoretical investigations. For example, the conventional synthesis route of polymeric carbon nitride, such as thermal polycondensation in air, typically yields a semicrystalline structure with mixed motifs, thereby complicating precise structural characterization.¹⁰ The reduced crystallinity of g-C₃N₄ is commonly attributed to incomplete polymerization and the persistence of hydrogen bonding within the synthesized sample, which results in a partially condensed framework of a melon-type g-C₃N₄ structure.¹¹

Recently, Negro et al.¹⁰ conducted a combined experimental and theoretical investigation on g-C₃N₄ structure using different characterization techniques such as XRD, FTIR, SEM, UV-Vis, and data from density functional theory (DFT) tools. Their study demonstrated that thermal polymerization of melamine produces a mixture of structural motifs comprising less-condensed melon frameworks together with highly condensed g-C₃N₄ domains. Another study by Zhang et al.¹² provides insights into the structural transformation mechanism of g-C₃N₄ synthesized via a molten salt-assisted procedure. They report a sequential phase evolution in which melon converts first into poly(heptazine imide) as an intermediate phase, and subsequently into poly(triazine imide) as the more stable phase. Collectively, earlier works have demonstrated that the structure, and consequently the properties, of g-C₃N₄ can be modulated via different synthesis routes such as precursor composition, condensation temperature, reaction rates and related factors.^9,13 Nonetheless, the fundamental question of what is the ground state structure of fully polymerized heptazine-based g-C₃N₄ has not been answered through experimental investigation.

In an effort to resolve the true ground state of pure heptazine-based g-C₃N₄, numerous theoretical investigations based on DFT and molecular dynamics have been conducted.^14–16 While numerous studies employ a graphene-like planar model of P6̄m2 symmetry for the g-C₃N₄ layer structure, other studies have allowed the geometry to flex into nonplanar structures. Recent theoretical work^11,17,18 has demonstrated that out-of-plane structures, such as buckled or corrugated layers, are energetically more favourable and offer a more realistic model of the fully polymerized heptazine phase. Moreover, Zhao et al.¹⁹ proposed three corrugated structures of P321, P3m1 and Pca2₁ symmetry for experimentally sythesized g-C₃N₄. These were found to be energicatlly more stable than the conventional flat P6̄m2 model. Other theoretical investigations that adopted a flat g-C₃N₄ model as the initial structure for adsorption or dopant studies consistently yielded buckled configurations, with planarity not being preserved once the symmetry is broken in the calculation.^20–23 These observations collectively suggest that out-of-plane buckling is an intrinsic feature of g-C₃N₄. However, it remains challenging, experimentally, to detect the subtly buckled heptazine phase. When the basal plane of the layered structure is viewed from above (with atomic-resolution microscopy techniques) the corrugated and planar g-C₃N₄ structures would be virtually indistinguishable.

Another important consideration in the search for the true ground state of heptazine-based g-C₃N₄ requires careful evaluation of how the two-dimensional (2D) layers are arranged to form a stacked three-dimensional (3D) bulk structure. The different stacking configurations between adjacent g-C₃N₄ sheets play a decisive role in determining the structural stability and electronic properties of the bulk phase. Although many theoretical studies^24–26 have modelled exfoliated g-C₃N₄ as an isolated monolayer to rationalize enhanced photocatalytic activity, the persistence of the (0002) diffraction peak indicates evidence of residual interlayer stacking. As a result, interlayer interactions and π–π overlap between adjacent sheets must not be neglected in describing the ideal phase of heptazine-based g-C₃N₄. In 2018, Sun et al.²⁷ investigated the electronic properties of graphitic carbon nitride with different building blocks and staggered between sheets in the bulk phase. In their study, triazine- and heptazine-based structures were examined under different stacking configurations (P6̄m2, R3m, Cmcm, and P6₃/mmc), and the results showed that bulk g-C₃N₄ polymorphs exhibit narrower band gaps than the corresponding monolayers, with the reduction attributed to interlayer band overlap between adjacent CN sheets. Other previous studies have also indicated that polymeric carbon nitride comprises stacked layers that preferentially adopt AB stacking registries as compared to AA.^18,28 However, the precise stacking order and local structural motifs remain uncertain due to experimental constraints, motivating further investigation.

A more recent study²⁹ demonstrated that the electronic and optical responses of polymeric g-C₃N₄ are dependent on the microstructural features, particularly the degree of condensation, out-of-plane corrugation, and interlayer interactions. These conclusions are consistent with prior work,^10,11 providing a strong rationale to systematically investigate the effect of buckling and stacking configurations in developing a more robust theoretical model for the fully condensed heptazine phase of g-C₃N₄.

Mridha et al.³⁰ recently advanced the earlier work of Im and colleagues²⁹ by providing a dynamic analysis of the ground-state structure of g-C₃N₄, systematically examining heptazine-, h-triazine-, and o-triazine-based polymorphs in both bulk and monolayer models. Their phonon calculations revealed that all planar polymorphs are dynamically unstable and relax to buckled structures, with heptazine identified as the energetically preferred structure.

While earlier studies have confirmed that buckling is intrinsic to the stability of g-C₃N₄, their scope remained limited to idealized monolayers and bulk phases, without assessing the influence of stacking registries on the relative stability of corrugated heptazine structures. To address this gap, the present study systematically investigates the interplay between buckling and different stacking configurations in heptazine-based g-C₃N₄, with a direct comparison between planar/buckled monolayers and the energetically preferred corrugated bulk phase. In addition, we have investigated how altering the assumed buckling and stacking structure affects the properties of phosphorus and nickel doping of g-C₃N₄. Our results reveal that the combination of buckling and appropriate stacking configurations stabilizes heptazine by yielding lower-energy corrugated structures.

2 Computational details

All periodic first-principles calculations performed in this study were carried out using plane-wave density functional theory (DFT) as implemented in the Vienna ab initio Simulation Package (VASP) code.^31,32 The interaction between ionic cores and valence electrons was described using the projector-augmented-wave (PAW) potentials provided with VASP.

Exchange–correlation energies were treated within the generalized gradient approximation (GGA) using the Perdew–Burke–Ernzerhof (PBE) functional. To properly account for weak interlayer interactions in the bulk models, van der Waals corrections were included using the DFT-D3 scheme³³ with Becke–Johnson damping. The PBE functional without dispersion corrections as well as the r²SCAN meta-GGA and HSE06 screened hybrid functionals were used for specific comparisons, as noted in the results presented below. These functionals were applied at the geometries that minimize the PBE-D3(BJ) energy. To assess the robustness of the predicted structure, one representative bulk structure was also fully optimized using different functionals (see Table S3 and Fig. S3). The plane wave energy cutoff was set as 500 eV, while the self-consistency energy convergence criterion was set to 10⁻⁶ eV. All structures were based on layers comprised of 2 × 2 supercells of the heptazine monomeric unit and an accompanying pore. The Brillouin zone was sampled using Γ-centered k-point meshes, specifically a 3 × 3 × 1 grid for 2D structures (incorporating a nominally 20 Å gap between layers) and a 3 × 3 × 7 grid for the 3D bulk systems. These grid point meshes were selected based on convergence tests that confirmed they are well-converged with respect to total energy (see Fig. S1).

To systematically search for lower-energy configurations, we introduced out-of-plane shifts of 0.5 Å for different atoms within heptazine to break the planar symmetry as depicted in Fig. 1. The monolayer was partitioned into four symmetry-related regions labelled a, b, c, and d. Each label corresponds to two pyridinic nitrogen atoms (N₂C sites), which are two-coordinated nitrogen atoms bonded to two carbon atoms within each heptazine unit. These N₂C atoms serve as the primary sites where out-of-plane perturbations were applied to initiate buckling. Model 1 represents the planar reference structure without any atomic perturbation, whereas models 2 to 5 incorporate systematic clockwise (+) and anti-clockwise (−) out-of-plane atomic shifts.

Fig. 1 Exploration of low-energy g-C₃N₄ buckling patterns initiated by atomic shifts/distortions of the planar monolayer.

To determine the stability of P and Ni dopants in the 2D/3D g-C₃N₄ structures, we carried out spin-polarized DFT calculations. The energy required to create a substitutional defect (E^sub_form) and an interstitial defect (E^int_form) were calculated using the following formation energy equations:

E^sub_form = E_doped − E_undoped − μ_dopant + μ_C/N (1)

E^int_form = E_doped − E_undoped − μ_dopant (2)

where E_doped and E_undoped are the total energies of the doped and pristine supercells, μ_dopant is the chemical potential of the added dopant atom, while μ_C/N corresponds to the chemical potential of the atom that occupied the substitution site prior to doping (either C or N). The chemical potential terms were taken as the per atom energy from the relevant pure phase: black phosphorous for P, metallic nickel for Ni, graphite for C and molecular dinitrogen (isolated in a large periodic box) for N. Care was taken to converge the k-point grids for these pure phase calculations (Fig. S7) to ensure that the corresponding uncertainties in eqn (1) and (2) due to sampling were significantly smaller than the values quoted in this work. A single dopant atom was used in each 2 × 2 supercell for the single layer stuctures and one dopant atom in the two layer supercell for the bulk structures.

The interlayer binding energy was calculated using the expression:^34,35

where A is the basal-plane area of the bulk unit cell, E_single-layer is the total energy of the optimized monolayer, and E_bulk is the total energy of the corresponding two-layer bulk supercell.

A limited investigation of molecular heptazine dimers was conducted using the Gaussian 16 code.³⁶ Two heptazine molecules (C₆N₇H₃) were brought together edge to edge in a constrained relaxation, controlling the N–N distance of one intermolecular pair and maintaining the planes of the two molecules at right angles to emphasise the N–N lone pair interaction. The MP2/def2-SVP level of theory was used for dimer calculations. Moreover, limited reverse Monte Carlo (RMC) simulations were performed on X-ray total scattering data from synthesised samples of g-C₃N₄.³⁷ These simulations were performed with the RMCProfile code.^38–40 Starting from a 10 000 atom supercell of planar bulk g-C₃N₄, the RMC iteratively perturbs the positions of these atoms until the calculated scattering pattern matches the experimentally measured data.

3 Results and discussion

3.1 Structural stability of g-C₃N₄ monolayer

We begin our discussion with the 2 × 2 supercell of the flat heptazine monolayer, constructed from its unit cell (a = b = 7.1 Å). Upon geometry relaxation using the PBE functional, the planar configuration persisted as an in-plane heptazine network with supercell lattice constants of a = b = 14.2 Å (Fig. 2a), consistent with previous theoretical work.^41,42 However, the optimized flat model is an unstable stationary point. Displacements in the positive and negative z directions are symmetry equivalent so that the z components of calculated forces remain exactly zero even though out-of-plane distortion lowers the energy. A series of controlled out-of-plane perturbations were therefore introduced to the planar geometry as shown in Fig. 1 to systematically search for lower-energy configurations. Each of these perturbed structures was subsequently relaxed to yield buckled configurations lower in energy than the planar reference (model 1), as summarized in Table S1. Among the corrugated structures considered, model 4 (defined by alternating clockwise and anti-clockwise atomic shifts in the pattern indicated in Fig. 1) emerged as the most stable, with the supercell about 2.9 eV lower in energy than the flat configuration at the PBE-D3(BJ) level. This energetic preference was consistent across other tested exchange–correlation functionals (r²SCAN and HSE06), providing strong evidence that the corrugated model offers the most realistic description of the ground state of single-layer heptazine-based g-C₃N₄. As shown in Fig. 2b, the optimized buckled configuration (model 4) was adopted as the representative corrugated model for heptazine monolayer in our study. As compared to the planar sheet (Fig. 2a), the buckled structure relaxed to a reduced lattice constant of 13.9 Å with a peak-to-peak buckling height of 1.35 Å, consistent with other studies.^16,30,43 Analysis of the structure after relaxation revealed that the planar model retained elongated bridge N–C bonds (∼1.47 Å), whereas the buckled structure allowed these bonds to contract slightly to ∼1.42 Å. The remaining N–C bonds increased slightly from 1.33 to ∼1.35 Å, reflecting redistribution of local strain.

Fig. 2 Relaxed 2 × 2 supercell of (a) planar and (b) buckled heptazine-based g-C₃N₄ monolayers.

Note that the energy difference between our model 4 and the flat model 1 is significantly larger than the differences found between the structures investigated by Mridha et al.³⁰ The periodicity of the single heptazine cell used by this earlier work limits the structures that can be found by their phonon-following approach to be of our model 2 type (cf. Fig. 1).

The shortening of the bridging C–N bonds upon buckling suggests buckling relieves steric hinderance between lone pairs of non-bonded nitrogens on the heptazine edges. This was investigated briefly by examining the potential energy profile of two heptazine molecules (C₆N₇H₃) brought together edge to edge in a constrained relaxation, controlling the nearest intermolecular N–N distance. Increasing the N–N distance from the 2.4 Å found in the flat structure to the 2.7 Å of the buckled geometry decreased the total energy by 0.43 eV at the MP2/def2-SVP level. As there are twelve such “intermolecular” N–N pairs in heptazine g-C₃N₄, this represents a significant reduction in energy through relief of these clashing electron-rich atoms. Although there will be some contribution from different levels of theory [MP2 versus DFT with PBE-D3(BJ)], that the total reduction of energy upon buckling is significantly less than 12 × 0.43 eV suggests that strain within the structure on buckling limits the distortion.

Our model 4 structure exhibits antiphase distortions (i.e. alternating distortion directions) along the 〈101̄0〉 and 〈011̄0〉 directions. This suggests energetic coupling between distortions similar to antiferroelectric coupling in ionic crystals. Such an interpretation would make the distortions frustrated on the hexagonal lattice of g-C₃N₄, as the symmetry-equivalent 〈1̄100〉 direction cannot simultaneously exhibit antiphase ordering in this finite 2 × 2 instantiation of the structure. Thus one may expect significant dynamical fluctuations of the out-of-plane distortion structure in real g-C₃N₄ crystals.

Our findings further support earlier reports^44,45 that corrugated structures resemble the flat g-C₃N₄ model (Fig. 2) when viewed from above the basal plane, which led to the widely used flat (metastable) P6̄m2 model in previous theoretical studies to interpret experimental data. However, having established the buckled configuration as the true energy minimum for the single-layer heptazine model, we proceed to examine the interlayer interactions and stacking registries of these corrugated sheets. The motivation to extend this structural investigation arises from experimental observations showing that exfoliated g-C₃N₄ nanosheets often retain interlayer spacings corresponding to the (0002) reflection in X-ray diffraction patterns. Moreover, while the precise interlayer registry of heptazine-based g-C₃N₄ remains unresolved, the stacking configuration cannot be neglected when accurately describing the synthesized heptazine sample. Therefore, in the following section, we systematically explore different stacking possibilities of heptazine-based g-C₃N₄ using first-principles theoretical calculations.

3.2 Stacking motifs in bulk g-C₃N₄

Similar to graphite, g-C₃N₄ consists of 2-dimensional layers held together by weak van der Waals forces, allowing adjacent sheets to take on a range of lateral displacements and orientations. This results in different stacking configurations, which constitute a fundamental structural variable that can be exploited to modulate the electronic and optical properties of g-C₃N₄.^46,47 We systematically investigated different stacking configurations of layered heptazine-based g-C₃N₄.

A series of initial configurations was tested to search for the most favourable interlayer configurations. In all cases we use a supercell comprising two layers of g-C₃N₄ sheets, yielding 2 × 2 supercells of the heptazine monomer. These are illustrated in Fig. 3.

Fig. 3 Structural investigation of different stacking motifs of heptazine-based g-C₃N₄ in both flat and buckled models to identify the most stable configuration. The asterisk (*) and subscripts c and f refer to 180° rotation, corrugated, and flat models, respectively. The side views are along the [112̄0] direction.

The bulk material with no lateral shift between each layer is denoted AA stacking. A range of AB-type stacking patterns were generated for testing. Those we label AB₁, AB₂ and AB₃ were generated by introducing an initial lateral shift of 0.6 Å along the a, −b or −2b − a directions, respectively, from the AA structure for every other layer. Another set of 180°-rotated AB variants was modelled, in which one flat or corrugated layer was first rotated by 180° relative to the perfectly aligned AA stacking. In this arrangement, the carbon atoms lie within the void regions, while only some nitrogen atoms, including the bridging nitrogen atoms (not carbon), overlap. As seen in Fig. 3, AB* is used to denote the structure resulting solely from this rotation, a structure that is usually denoted simply AB in other work. We also consider applying similar lateral shifts to those described above from this AB* structure, denoted , and . This notation was extended with an additional subscript f or c to distinguish between structures constructed from flat or corrugated g-C₃N₄ layers, respectively. For example, is a staggered model comprising corrugated layers of alternating orientation that has been initially shifted along the b axis.‡ All structures were relaxed with the PBE-D3(BJ) functional.

The optimized lattice parameters, lattice angles, and relative energies [computed using various exchange–correlation (XC) functionals] for each configuration are presented in Table 1. Consistent with previous theoretical studies,^27,28 the flat AA_f reference model (with lattice constants a = 14.2 Å and c = 6.73 Å) was the highest energy structure. However, introducing a 180° rotation in the flat AB_f stacking motif, denoted as , lowers the energy by ∼3.1–4.0 eV across all the XC functionals, highlighting that changes in the stacking registry can significantly influence structural stability. In contrast, the corrugated AA_c structure undergoes an in-plane contraction (a = 13.9 Å) and a slight expansion along c (∼ 7.18 Å) and lowers the energy by ∼6.6 eV. The enhanced stability of the AA_c stacking configuration, relative to the previously discussed models, can be attributed to the relief of interlayer electrostatic repulsion through out-of-plane buckling.

Table 1 Calculated lattice parameters [PBE-D3(BJ)] and relative stabilities (ΔE, eV per supercell) of heptazine stacking motifs from different XC-functionals

Stacking motifs	Lattice constants (Å)	Lattice angles	Relative energy, ΔE (eV)
			PBE	PBE-D3(BJ)	r^2SCAN	HSE06
AAf	a = b = 14.2, c = 6.73	α = β = 90°, γ = 120°	0.00	0.00	0.00	0.00
	a = b = 14.2, c = 6.38	α = β = 90°, γ = 120°	−3.13	−4.02	−3.81	−3.44
AAc	a = b = 13.9, c = 7.18	α = β = 90°, γ = 119.8°	−6.61	−6.59	−6.91	−6.62
AB1c	a = b = 13.9, c = 6.97	α = β = 89.8°, γ = 119.8°	−7.68	−7.89	−8.16	−7.79
AB2c	a = b = 13.9, c = 6.96	α = β = 90°, γ = 119.9°	−7.69	−7.89	−8.15	−7.79
AB3c	a = b = 13.9, c = 6.97	α = β = 90°, γ = 119.9°	−7.70	−7.90	−8.18	−7.81
	a = b = 14.0, c = 7.60	α = 92°, β = 87.9°, γ = 120.1°	−7.32	−6.14	−6.52	−7.20
	a = b = 14.0, c = 6.90	α = β = 90°, γ = 120.2°	−6.73	−6.93	−7.16	−6.83
	a = b = 14.0, c = 6.91	α = β = 89.8°, γ = 119.8°	−6.73	−6.92	−7.14	−6.82
	a = b = 13.5, c = 6.98	α = β = 89.8°, γ = 119.8°	−8.87	−9.64	−9.96	−9.29

---

Modifying the interlayer arrangement in the corrugated bulk structure results in a further reduction in total energy, with the AB_3c registry attaining the lowest relative energy (−7.81 eV at the HSE level) within this stacking series. Moreover, among the 180°-rotated corrugated stacking series, shows moderate stability, whereas the addition of translational shifts in achieves the lowest energy stacking configuration across all the XC functionals. Consistent with these structural relaxations, the corrugated stacking registries exhibit buckling heights ranging from ∼0.97 to 2.84 Å, including 1.28 Å for AAc and 0.97 Å for , while the AB_1c–AB_3c motifs show distortions of ∼1.31–1.42 Å, reaching up to 2.84 Å for .

To further investigate the interlayer coupling strength of the various stacking registries in g-C₃N₄ the interlayer binding energy is considered. This provides a quantitative measure of the cohesion between adjacent heptazine sheets.

The interlayer binding energy is calculated using eqn (3). Table S2 provides details of the interlayer binding energies for all planar and corrugated stacking registries. Across these configurations, the binding energy values fall within 6–17 meV Å⁻², lower than the binding/exfoliation energies of some other van der Waals materials such as graphene (22.5 meV Å⁻²)⁴⁸ and MoS₂ (26.3 meV Å⁻²).³⁴ This shows that all g-C₃N₄ stacking registries considered here are weakly bound and readily exfoliable, yet the relative differences among them remain significant. For example, the AA_f registry exhibits the lowest interlayer binding, while the configuration shows the highest binding energy, indicating a more favorable alignment that enhances structural stability. This trend highlights that even minor registry shifts or corrugation patterns can modulate the weak van der Waals interactions between adjacent heptazine layers.

Overall, the results of our theoretical investigation show the profound influence of the synergistic effect of corrugation and stacking registry on determining the true ground state of layered heptazine-based g-C₃N₄. The emergence of these newly identified corrugated stacking configurations (all exhibiting P1 symmetry) provides a robust structural framework for systematically investigating how stacking transitions through, for example, applying mechanical shear, can modulate the electronic properties of bulk g-C₃N₄, as will be discussed in the next section.

3.3 Electronic structure of pristine heptazine models

To gain insight into how structural corrugation modulates the electronic properties of heptazine-based g-C₃N₄, we investigated the electronic band structures and partial densities of states (PDOS) of both the planar and buckled monolayers using PBE and hybrid HSE06 functionals. As depicted in Fig. 4, the computed band gaps, E_g, characterized by the energy difference between the valence band maximum (VBM) and conduction band minimum (CBM), for the flat monolayer are 1.25 eV (PBE) and 2.82 eV (HSE06), consistent with previously reported literature values.^16,49–51 In contrast, buckling increases the band gap of the corrugated monolayer to 1.87 eV (PBE) and 2.97 eV (HSE06). The out-of-plane displacement in the buckled structure induces angular distortions in the C–N bonds, which lower the local symmetry and reduce orbital overlap between adjacent units, thereby leading to a wider band gap. As expected, the PBE functional underestimates the band gap due to the many-electron self-interaction error, but the HSE06 hybrid functional provides a band gap closer to expectations from experiment. However, the overall electronic trend of the band gap increasing from the planar to the buckled configuration remains consistent across both exchange–correlation functionals, confirming that the intrinsic structural buckling is a fundamental feature that influences the electronic structure of g-C₃N₄.

Fig. 4 Electronic band structures and partial density of states (PDOS) of the planar (left) and buckled (right) heptazine-based g-C₃N₄ monolayers. Upper panels PBE, lower panels HSE06. Energies are shown relative to the Fermi level reported by Vasp.

Fig. 4 also shows the corresponding PDOS of both planar and buckled configurations, revealing the contributions of each atomic orbital type to the electronic band. In the valence band (VB), the major contribution originates mainly from the N-2p lone-pair orbital, while the conduction band (CB) is composed of the hybridization of C-2p and N-2p states. Moreover, both the band structure and PDOS of the buckled monolayer reveal an upward shift of the conduction-band edge compared to the planar model, resulting in a larger band gap (although we note that the Fermi energy is not defined meaningfully in calculations of this type). Our observations are in good agreement with previous reports.^45,52,53

We now shift our attention to examining how different stacking motifs of the heptazine units influence the electronic properties of bulk g-C₃N₄. The layered structure of 3D–g-C₃N₄ introduces additional band dispersion along the out-of-plane direction, resulting in a narrower band gap compared to the isolated monolayer. For example, comparing the band structure plots in Fig. 4 and 5, all the bulk stacking configurations exhibit smaller band gaps than the monolayer. For the flat models (at the PBE level), the band gap decreases from 1.25 eV in the monolayer to 0.59 eV (AA_f) and 1.03 eV , and from 2.82 eV to 1.93 eV (AA_f) and 2.06 eV at the HSE06 level. A similar trend is also observed for buckled models, where the corresponding band gaps reduce from 1.87 eV in the monolayer to 1.59 eV (AA_c) and 1.60 eV at the PBE level, and from 2.97 eV to 2.66 eV (AA_c) and 2.65 eV using the HSE06 hybrid functional.

Fig. 5 Electronic band structures of selected heptazine-based bulk g-C₃N₄ obtained using PBE and HSE06 are shown, with correspoding top and side views of the structures included. Each panel corresponds to a different stacking motif: (a) AA_f, (b) AB^*_f, (c) AA_c, and (d) AB^*_c, where subscripts f and c denote flat and corrugated models, respectively.

Next, we extend our analysis to compare how different stacking arrangements (in the flat and corrugated models) influence the electronic response of bulk g-C₃N₄. It can be seen from Fig. 5a–d and S2 that introducing a 180° rotation in the AA stacking registry, together with the presence of buckling in either stacking configuration, modulates the bulk electronic properties. While rotating the flat AA stacking motif moderately widens the band gap from 1.93 eV to 2.06 eV (HSE06), the introduction of buckling in the AA stacking (denoted as AA_c), even without rotation, results in a wider band gap of 2.66 eV (HSE06), which is in close agreement with experimental observations (2.7 eV).⁵⁴ Moreover, applying a 180° rotation to the buckled AA stacking induces minimal changes in the band gap, indicating that structural buckling exerts a more dominant role in determining interlayer electronic coupling compared with rotational misalignment alone. These findings suggest that lateral translation or combined rotation-translation of the layers, for example by applying shear stress to the material, could serve as a viable strategy to tune the bulk electronic properties of van der Waals layered g-C₃N₄. As shown in Fig. S2, the electronic response of the different corrugated stackings modulates the band dispersion and PDOS of bulk g-C₃N₄.

A more detailed comparative analysis of the band gaps in Fig. 6 quantifies the effect of these corrugated stacking variants. Relative to the unshifted corrugated models (i.e. AA_c, ), lateral translation increases the band gap by about 6%. In contrast, the combined rotation–translation raises it by roughly 15%, with the most stable structure, , exhibiting the widest band gap (3.05 eV) in these stacking series. Our findings reveal that while stacking of flat layers reduces the band gap, corrugation of stacked layers both enhances the thermodynamic stability and yields larger band gaps, consistent with the postulation of Im et al.¹¹ Accordingly, one may speculate that the band gap modulation approach demonstrated in this work may be extended to other two-dimensional materials to control their electronic properties, possibly including charge transport and photoresponse characteristics.

Fig. 6 Relative band-edge positions and band gap variation (HSE06) across different stacking configurations of heptazine-based g-C₃N₄. The Fermi level, E_f, reported by Vasp is set to zero. The asterisk (*) and subscripts ‘f’ and ‘c’ refer to 180° rotation, flat and corrugated models, respectively.

3.4 Dopant roles in monolayer and bulk g-C₃N₄

To investigate the influence of a dopant on 2D planar and 2D/3D buckled g-C₃N₄—and more importantly, how this is affected by the variations in structure discussed above—we select representative non-metal (P) and metal (Ni) dopants. Starting with the 2D buckled monolayer, we examined all possible substitutional P sites, at C1–C6 and N1–N7 (Fig. 2), as well as the Ni interstitial site, using the defect formation energy expressions defined in eqn (1) and (2). Our calculations reveal that among the C dopant sites, P@C6 is the most thermodynamically favourable configuration, with a formation energy of 0.19 eV (see Fig. 7 and S5). The P@N sites require significantly higher energies (0.60–1.90 eV, Fig. S6).

Fig. 7 Structural configurations and formation energies of P and Ni dopants in planar and corrugated monolayer (2D) and bulk (3D) g-C₃N₄. (a–c) P-doped and (d–f) Ni-doped systems, corresponding to 2D planar, 2D corrugated, and 3D corrugated structures, respectively.

The P@C6 site was used as a model defect to examine how dopant accommodation varies with buckling and stacking structure, comparing the 2D-planar and buckled monolayers with the buckled 3D-stacked g-C₃N₄. Fig. 7 shows that the formation energy changes dramatically from the planar sheet (1.64 eV) to the buckled monolayer (0.19 eV), whereas the 3D-buckled structure exhibits an intermediate stability (0.67 eV). The variation in formation energies, despite using the same dopant sites (P@C6) across all structures, indicates that each host lattice accommodates the dopant differently. In the planar monolayer (Fig. 7a), the P dopant is constrained by the in-plane symmetry of the lattice, resulting in a P–N bond length of 1.59 Å, longer than the intrinsic C–N bond length of 1.33 Å, giving rise to local strain. In contrast, the buckled monolayer alleviates this strain by allowing the dopant to relax out of plane, as reflected by the P–N bond length (1.63–1.64 Å, Fig. 7b). In the 3D-buckled structure (Fig. 7c), the local bonding geometry around the P dopant is nearly identical to that of the buckled monolayer, with comparable P–N distances. However, the difference in the formation energies between the 2D buckled monolayer and the 3D buckled lattice can be attributed to the long-range interlayer interactions in the bulk and not from the local bonding environment. Similarly for the Ni dopant at the same interstitial sites across the planar, buckled monolayer, and 3D-buckled structures (Fig. 7d–f), the formation energies vary with 1.63 eV for planar 2D, 1.72 eV for buckled 2D, and 1.35 eV for buckled 3D. The variation in dopant energetics between P and Ni at identical sites across flat, buckled, and layered g-C₃N₄ structures demonstrates that accurate theoretical modelling of g-C₃N₄ must account for intrinsic structural stability as well as the influence of interlayer coupling.

The buckling and stacking structure influenced the effect of dopants on the calculated electronic structure. Fig. 8a–c and S8a–c show that doping with P does not significantly alter the bulk band structure of the underlying g-C₃N₄ structure. The significant change to the electronic structure is the introduction of mid-gap states associated with the P dopant atom. For the neutral systems modelled here, a single electron occupies one of these states, with the opposite spin state remaining unoccupied. Clearly, such mid-gap states could be expected to have a significant effect on the photocatalytic properties of doped g-C₃N₄.

Fig. 8 Spin-polarized projected density of states (PDOS) of P- and Ni-doped g-C₃N₄; (a–c) P-doped and (d–f) Ni-doped systems across 2D planar, 2D buckled, and 3D buckled configurations.

The position of these states within the gap, and the separation between the occupied and unoccupied states, are strongly affected by the modelled structure. The planar 2D structure yields a small separation between the occupied and unoccupied states, with both lying close to the conduction band of g-C₃N₄. Allowing the monolayer to buckle out of the plane separates the occupied and unoccupied states by around 1 eV, with the occupied state lying close to the valence band. Stacking g-C₃N₄ into a bulk structure causes both mid gap states to move higher in the band gap.

Charge density difference (CDD) and Bader charge transfer calculations for these P-doped structures are represented in Fig. 9a–c. These again demonstrate significant differences between the different g-C₃N₄ parent structures, with significant differences between the amount of charge transferred to the g-C₃N₄ structure depending on whether the doped parent is planar, buckled, or in bulk.

Fig. 9 Charge density difference (CDD) distributions and Bader charge variations (Δq) of P- and Ni-doped g-C₃N₄ in planar 2D, buckled 2D, and buckled 3D structures: (a–c) P-doped and (d–f) Ni-doped systems. The yellow and cyan regions denote electron accumulation and depletion, respectively. The isosurface is taken as 0.003 e bohr⁻³.

For g-C₃N₄ doped with interstitial Ni, the parent band structure was again largely unchanged from the undoped material's band structure for the monolayer cases (Fig. 8d, e and S8d, e). Occupied Ni 4p states were introduced at the top of the gap for the flat structure, with Ni 3d lying at the bottom of the gap. The wider gap of the more stable buckled monolayer again presented occupied Ni 4p states at the top of the gap, but exposed Ni 3d states lying in the parent material band gap at slightly destabilized energies. Introducing a Ni dopant to the bulk AB_3c structure (Fig. 8f and S8f) opened the parent material band gap to almost 2 eV, with occupied Ni 3d states again appearing in the lower half of the parent gap. In this case, occupied Ni 4p states were stabalised to similar energies to the Ni 3d states, maintaining a nominal intrinsic band gap of around 0.7 eV. Fig. 9d–f indicates that the electron transfer from the Ni dopant was similar in the planar 2D and buckled bulk cases but reduced by more than 20% when the monolayer was allowed to buckle to a more stable configuration.

3.5 Experimental verification

To this point, this paper, along with others, presents a clear case based on DFT that the “graphitic” layers of g-C₃N₄ are not in fact flat planes like graphite, but rather strongly buckled structures. Although modern DFT approaches are reliable, sceptics may question whether this behaviour is an artefact of the sophistication of the exchange–correlation functional, the treatment of dispersion, or other deficiencies of the theory.

As a final independent piece of evidence that the layers of g-C₃N₄ are indeed buckled away from planarity, we present some experimental evidence. X-ray total scattering data for g-C₃N₄ was simulated via reverse Monte Carlo calculation with the RMCProfile code.^38,39 During the RMC run, the atoms of g-C₃N₄ moved consistently away from planar structures in order to give a good fit to the X-ray scattering data, which could not be obtained from a structure comprised of flat planes of atoms. A side view of this real-space structure is shown in Fig. 10.

Fig. 10 (a) and (b) Experimental X-ray total scattering G(r) (black) compared with planar and buckled AB-stacked g-C₃N₄ models (red), showing poor agreement for the planar structure (R_w = 51.49%) and markedly improved agreement for the buckled AB model (R_w = 15.80%); dashed curves denote C–C, C–N, and N–N partial correlations. (c) Ball-and-stick representation of the RMC-refined g-C₃N₄ structure fitted to X-ray total scattering data, showing part of the 10 000-atom supercell viewed along [112̄0]. (d) Comparison of the experimental G(r) (black) with the final RMC model (red).

The out-of-plane displacements generated by the RMC simulation are similar in magnitude to those exhibited by our DFT structures. Note that the RMC structures maintained atomic positions consistent with the structure as a heptazine polymer throughout. Further details of these and related RMC simulations will be presented in a related paper.³⁷

4 Conclusions

In this study, we systematically explore the structural stability of heptazine-based g-C₃N₄ from a monolayer to different stacking configurations and how the presence of dopants in these configurations modulates the structural and electronic properties. Our results show that:

• The widely assumed planar (P6̄m2) model of g-C₃N₄ is neither energetically stable nor a reasonable approximation to the structure of ordered g-C₃N₄. Allowing out-of-plane atomic relaxation yields a lower-energy buckled structure that alleviates the in-plane N-lone-pair repulsion inherent to the planar model, whereas doping with heteroatoms is not required to induce buckling.

• In the layered heptazine-based g-C₃N₄ structure, the stacking patterns previously demonstrated to be the lowest in energy for flat layers are significantly higher in energy than other stacking patterns when the layers are buckled.

• The electronic structure of g-C₃N₄ is sensitive to both buckling and stacking registry. Although buckling generally widens the band gap relative to the planar monolayer, different corrugated stacking configurations yield distinct band gap values, which could potentially be an effective mechanism for tuning the electronic properties of bulk g-C₃N₄. For example, the corrugated stackings (P1 symmetry) uncovered in this study display varied electronic responses, with the band gaps of certain stacking configurations aligning closely with experimentally reported values. The most stable structure found in this work was the one labelled while AB_3c showed the best match to the experimental band gap (calculated at the HSE06 level).

• The calculated stability and electronic effect of representative metal and non-metal single atom dopants is strongly impacted by the treatment of both buckling and stacking.

We present a number of corrugation patterns that were calculated to be lower in energy than any of those previously reported. Key to these lower energy structures was the use of a 2 × 2 supercell within each plane, allowing adjacent heptazine units to buckle in opposite directions. The hexagonal structure of g-C₃N₄ suggests this alternation of structure may be frustrated in the third equivalent in-plane direction.

Within these corrugated structures, P and Ni dopants exhibit markedly different electronic responses compared with the flat model, demonstrating that dopant functionality depends primarily on the host structural motif. As such, reliable property prediction and rational design of both pristine and doped g-C₃N₄ require careful consideration of accurate structural models to enable the development of high-performance g-C₃N₄ photocatalysts. In particular, studying the effect of dopants by using DFT calculations based on either flat or buckled g-C₃N₄ monolayers is likely to be misleading if the target material is a synthesisable (i.e. multilayer or bulk) doped g-C₃N₄. Such calculations yield incorrect predictions of both the propensity of dopants to enter the g-C₃N₄ lattice and the effects of such dopants on the electronic (and hence photocatalytic) properties of the doped material.

Conflicts of interest

There are no conflicts of interest to declare.

Data availability

The data supporting the findings of this study have been included in this article and in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6ma00065g.

Acknowledgements

This work was supported by the Australian Research Council Centre of Excellence for Carbon Science and Innovation (CE230100032), the Australian Laureate Fellowship program (FL210100017), and the Linkage Project (LP200200472). This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI Australia), an NCRIS-enabled capability supported by the Australian Government. Access to NCI resources was facilitated by the UNSW HPC Scheme (https://doi.org/10.26190/PMN5-7J50).

Notes and references

1. Y. Tong, J. Xia, Y. Hu, Y. He, G. He and H. Chen, Chem. Commun., 2025, 61, 1509–1532.
2. H. M. Tofil, R. Ghazi, F. A. Ghaib, R. N. Dara, I. Kebaili, I. Boukhris, H. Ding and Z. Rehman, Sustainable Energy Fuels, 2025, 9, 2900–2927.
3. O. K. Abass, G. F. Abass, Y. Lou, O. E. Ajayi, W. J. Onyia, P. H. Bassey, O. T. Faloye and D. Bala, Curr. Res. Green Sustainable Chem., 2025, 11, 100501.
4. W. Fang, D. Wu, W. Dai, S. Wang, B. Hu, G. Chen and S. Yu, Desalination, 2025, 614, 119194.
5. M. A. Khan, S. Mutahir, I. Shaheen, Y. Qunhui, M. Bououdina and M. Humayun, Coord. Chem. Rev., 2025, 522, 216227.
6. X. Zhao, Q. Liu, X. Li, H. Ji and Z. Shen, Chin. Chem. Lett., 2023, 34, 108306.
7. P. Phuong Ly, D. V. Nguyen, T. Anh Luu, M. Chien Nguyen, P. Duc Minh Phan, H. Phuoc Toan, T. Van Nguyen, M. T. Pham, T. Dieu Thi Ung, D. Danh Bich, H. Thi Pham, H. Thi Ngoc Nguyen, W. Jong Yu, S. Hyun Hur, N. Quang Hung and H. T. Vuong, Chem. Eng. J., 2025, 504, 158657.
8. M. Rizwan, A. Zada, S. Azizi, P. Fazil, M. Naeem, G. Murtaza, N. A. Ibrahim, N. S. Basher, Z. Liu and M. H. Eisa, Int. J. Hydrogen Energy, 2025, 97, 1126–1152.
9. K. Zhang, S. Tian, X. Wang, B. Li, J. Pang, J. Xu, D. Yang, F. Dai, Z. Wu, X. Chen, X. Wang, L. Wang and J. Xing, Green Chem., 2025, 27, 3863–3868.
10. P. Negro, F. Cesano, S. Casassa and D. Scarano, Materials, 2023, 16, 3644.
11. C. Im, B. Kirchhoff, I. Krivtsov, D. Mitoraj, R. Beranek and T. Jacob, Chem. Mater., 2023, 35, 1547–1559.
12. K. Zhang, C. Liu, Q. Liu, Z. Mo and D. Zhang, Catalysts, 2023, 13, 717.
13. D. Ma, Z. Zhang, Y. Zou, J. Chen and J. W. Shi, Coord. Chem. Rev., 2024, 500, 215489.
14. Y. Wang, C. Wu, Y. Tian, L. Yan, H. Tan and Z. Su, Appl. Surf. Sci., 2018, 453, 442–448.
15. J. Wang, D. Hao, J. Ye and N. Umezawa, Chem. Mater., 2017, 29, 2694–2707.
16. Q. Gao, X. Zhuang, S. Hu and Z. Hu, J. Phys. Chem. C, 2020, 124, 4644–4651.
17. G. D. Fao and J. C. Jiang, Mol. Catal., 2022, 526, 112402.
18. G. O. Hartley and N. Martsinovich, Faraday Discuss., 2021, 227, 341–358.
19. L. Zhao, X. Shi, J. Li, T. Ouyang, C. Zhang, C. Tang, C. He and J. Zhong, Phys. E, 2021, 128, 114535.
20. L. M. Azofra, D. R. MacFarlane and C. Sun, Phys. Chem. Chem. Phys., 2016, 18, 18507–18514.
21. X. Zhou, C. Zhao, C. Chen, J. Chen and X. Chen, Diamond Relat. Mater., 2022, 125, 108995.
22. W. Zhang, Z. Zhang, S. Kwon, F. Zhang, B. Stephen, K. K. Kim, R. Jung, S. Kwon, K. B. Chung and W. Yang, Appl. Catal., B, 2017, 206, 271–281.
23. Y. Ji, H. Dong, H. Lin, L. Zhang, T. Hou and Y. Li, RSC Adv., 2016, 6, 52377–52383.
24. C. Yang, Z.-Y. Zhao, H.-T. Wei, X.-Y. Deng and Q.-J. Liu, RSC Adv., 2021, 11, 4276–4285.
25. G. Fu, W. Zhen, H. Wang, X. Zhou, L. Yang and J. Zhang, Molecules, 2024, 29, 5339.
26. D. K. Gorai, S. K. Kuila, A. Kumar, M. I. Ahmad and T. K. Kundu, Appl. Surf. Sci., 2023, 623, 157031.
27. S. P. Sun, S. Gu, J. H. Sun, F. F. Xia and G. H. Chen, J. Alloys Compd., 2018, 735, 131–139.
28. S. Zuluaga, L.-H. Liu, N. Shafiq, S. M. Rupich, J.-F. Veyan, Y. J. Chabal and T. Thonhauser, Phys. Chem. Chem. Phys., 2015, 17, 957–962.
29. C. Im, M. M. Elnagar, B. Kirchhoff, D. Mitoraj, I. Krivtsov, A. Farkas, R. Beranek and T. Jacob, J. Mater. Chem. C, 2025, 13, 8682–8693.
30. M. M. Mridha, S. Kellici and J. Buckeridge, J. Phys. Chem. C, 2025, 129, 15109–15121.
31. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
32. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
33. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
34. Z. Bian, J. Miao, Y. Zhao and Y. Chai, Acc. Mater. Res., 2022, 3, 1220–1231.
35. J. H. Jung, C.-H. Park and J. Ihm, Nano Lett., 2018, 18, 2759–2765.
36. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, et al., Gaussian 16 Revision C.01, 2016.
37. A. Ahmed, G. Abass, M. T. Naing, C. Yang, D. Gunawan, H. Huang, Q. Zhai, Z. Zhang, C. Xie, T. Lu, N. Cox, L. Dai, R. Amal, T. Frankcombe and Y. Liu, Atomically coordinated nonmetal co-doped g-C₃N₄ for efficient H₂O₂ photosynthesis, Appl. Catal. B-Environ. Energy, 2026, submitted.
38. M. G. Tucker, D. A. Keen, M. T. Dove, A. L. Goodwin and Q. Hui, J. Phys.: Condens. Matter, 2007, 19, 335218.
39. Z. Liu, H. Peng, T. Lu, T. Wu, C. Yang, Z. Fu, Z. Hong, J. Xie, T. Honda, Y. Chen, W. Hu, F. Xu, Z. Lin, Y. Liu, S. Zhang, G. Wang and J. Chu, J. Am. Chem. Soc., 2025, 147, 41620–41628.
40. H. Wang, Y. Liu, L. Li, J. Zhang, M. Luo, Z. Sun, Y. Xiao, X. Zhai, L. Wu, H. Zhang, B. Ye, C. Yang, X. Zhang and M. Zhou, Angew. Chem., Int. Ed., 2025, 64, e202506960.
41. Q. Gao, S. Hu, Y. Du and Z. Hu, J. Mater. Chem. A, 2017, 5, 4827–4834.
42. H.-Z. Wu, L.-M. Liu and S.-J. Zhao, Phys. Chem. Chem. Phys., 2014, 16, 3299–3304.
43. X. Liu, W. Kang, W. Zeng, Y. Zhang, L. Qi, F. Ling, L. Fang, Q. Chen and M. Zhou, Appl. Surf. Sci., 2020, 499, 143994.
44. L. Zhao, X. Shi, J. Li, T. Ouyang, C. Zhang, C. Tang, C. He and J. Zhong, Phys. E, 2021, 128, 114535.
45. B. Mortazavi, M. Shahrokhi, F. Shojaei, T. Rabczuk and X. Zhuang, Comput. Mater. Today, 2025, 5, 100024.
46. C. Bao, W. Yao, E. Wang, C. Chen, J. Avila, M. C. Asensio and S. Zhou, Nano Lett., 2017, 17, 1564–1568.
47. C. Fox, Y. Mao, X. Zhang, Y. Wang and J. Xiao, Chem. Rev., 2024, 124, 1862–1898.
48. Y. Huang, Y.-H. Pan, R. Yang, L.-H. Bao, L. Meng, H.-L. Luo, Y.-Q. Cai, G.-D. Liu, W.-J. Zhao, Z. Zhou, L.-M. Wu, Z.-L. Zhu, M. Huang, L.-W. Liu, L. Liu, P. Cheng, K.-H. Wu, S.-B. Tian, C.-Z. Gu, Y.-G. Shi, Y.-F. Guo, Z. G. Cheng, J.-P. Hu, L. Zhao, G.-H. Yang, E. Sutter, P. Sutter, Y.-L. Wang, W. Ji, X.-J. Zhou and H.-J. Gao, Nat. Commun., 2020, 11, 2453.
49. R. Ahmed and A. K. Manna, J. Phys. Chem. C, 2022, 126, 18208–18215.
50. B. Mortazavi, F. Shojaei, M. Shahrokhi, M. Azizi, T. Rabczuk, A. V. Shapeev and X. Zhuang, Carbon, 2020, 167, 40–50.
51. B. Liang, Y. Rao and X. Duan, RSC Adv., 2019, 9, 38724–38729.
52. Q. He, F. Zhou, S. Zhan, Y. Yang, Y. Liu, Y. Tian and N. Huang, Appl. Surf. Sci., 2017, 391, 423–431.
53. B. Zhu, L. Zhang, B. Cheng and J. Yu, Appl. Catal., B, 2018, 224, 983–999.
54. L. Chen, M. A. Maigbay, M. Li and X. Qiu, Adv. Powder Mater., 2024, 3, 100150.

---

Footnotes

† These authors contributed equally to this work.

‡ Note that, at the time of writing, the structure labelled here as is listed as the stable structure in the Materials Project online database.

---