A practical computational protocol for photocatalytic reactions beyond ground-state approximations

Abstract

Theoretical studies of heterogeneous photocatalysts typically discuss the band level alignment obtained by ground-state DFT, which does not capture the physics of light-driven processes. Here, we present a new computational protocol in which excited states are explicitly considered in the Gibbs free energy diagrams. This is applied to prototypical H₂O oxidation and O₂ reduction reactions on a single-atom Co-embedded g-C₃N₄ cocatalyst.

---

Photocatalytic energy conversion with transition metal single-atom cocatalysts (TM-SACs) includes applications in photocatalytic H₂ evolution, CO₂ reduction, nitrogen fixation, and organic synthesis.¹ Recently, outstanding progress has been achieved in H₂O₂ production by O₂ reduction and H₂O oxidation.^2–8 In this exploratory study, we use the Co/g-C₃N₄ TM-SAC as a prototypical system to evaluate possible pathways within the oxygen evolution reaction (OER) network:

• standard 4e⁻ OER and 2e⁻ H₂O₂ formation through the classical *O and nonclassical *(OH)₂ intermediates,

• oxygen reduction reaction (ORR) to H₂O₂ (2e⁻) and to H₂O (4e⁻).

Reaction intermediates and pathways are depicted in Fig. 1a. Atomically dispersed cobalt in graphitic carbon nitride (g-C₃N₄) was experimentally found to improve H₂O₂ formation from the H₂O rate constant,² reduce the activation barriers of water oxidation and oxygen reduction,³ and enhance H₂O₂ production from O₂, despite highly positive reaction steps calculated from ground-state DFT.⁴ According to previous theoretical investigations of Co/g-C₃N₄, we placed Co in the plane of the heptazine-based 2D carbon nitride network in an M–N₄ bonding environment.^6,9

Fig. 1 (a) Schematic of the reaction pathway; orange and purple colors correspond to the 4e⁻ O₂ evolution and the 2e⁻ H₂O₂ evolution pathways, respectively. (b) Computational protocol hierarchy used to compute G* – excited free energy profiles. At the topmost Level 0, the ground-state picture is presented (E_EXC = 0, G* = G). Towards the bottom, increasingly robust and computationally intensive approximations are summarized with different quantities used as E_EXC: E_g (Kohn–Sham band gap), VEE (vertical excitation energy), and AEE (adiabatic excitation energy).

Here, we present our novel computational protocol (Fig. 1b), which allows us to find the correct theoretical approach needed to capture, at the DFT level, the preference towards photocatalytic H₂O₂ production in Co/g-C₃N₄.^2–4 In this work, we describe purely photocatalytic processes at room temperature, i.e. without an applied external potential. In this description, we use Kohn–Sham (KS) DFT to obtain the Gibbs free energy G₀ = E_tot + E_ZPE + H − TS = E_tot + G(T), where E_tot is the total energy of the optimized system, E_ZPE is the zero point energy, H is the enthalpic term, and TS is the entropic term. The reaction thermodynamic profile is expressed in Gibbs free energy difference (ΔG) from one step to the next. The overall ΔG values are presented in Table S1, with each contribution listed separately in Table S2. The reaction mechanisms consist of proton-coupled electron transfer (PCET) steps with free energy described by the relation . The reaction thermodynamics are evaluated in both ground and excited states. Since the scope of this study is to assess the catalytic activity of the chosen reaction using the proposed computational protocol, a few assumptions have been considered: (1) light absorption has the correct wavelength to excite the substrate to the conduction band; (2) the generated hole–electron pair is transferred at the interface; (3) the charges are fully separated and reach the catalytic center without recombination (verified by examining the charge density difference between the excited and ground states, see Fig. S3).

The standard ΔG profiles are expanded by explicitly considering the excited-state contributions in the Gibbs free energy calculations. In the proposed protocol, we consider an additional term E_EXC, in order to obtain the total expression for intermediate free energy G* = G₀ + E_EXC (here the ‘*’ denotes the excited state treatment and index ‘0’ denotes the ground state). The treatment of E_EXC considers different approximations for introducing the excited state contributions. In the ground state, E_EXC simply equals 0. We call this the ‘Level 0’ description. At Level 0+, ground-state calculations are still applied, but E_EXC = E_CBM − E_VBM = E_g, the Kohn–Sham eigenvalue band gap of the system. This represents a crude approximation for the energy increase when the material is illuminated by light. At Level 1, the vertical excitation energy E_EXC = VEE is considered, calculated using the ΔSCF method.¹⁰ This is obtained by comparing the total energies of excited and ground-state systems without relaxing the atomic geometry of the excited state. Finally, at Level 2, we allow excited systems to relax to their adiabatic excitation energy E_EXC = AEE. The ΔSCF procedure used at the two topmost levels involves fixing band occupations at a desired configuration (‘non-Aufbau’), instead of filling them from the lowest to the highest energy bands (the ‘Aufbau’ principle). The density optimization algorithm is chosen so that it is possible to proceed without diagonalizing the Kohn–Sham matrix, which would inherently reorder the orbitals by increasing single-particle energy. This ensures converging to the desired excited state, instead of the ground state (for exact details, see the SI). Here, the selected configuration involves depopulating the one-electron valence band maximum (VBM) and filling the conduction band minimum (CBM) of the g-C₃N₄ support material.

Each of these modifications to the Gibbs free energy state function G describes slightly different physics, which we will outline in more detail. At the Level 0 approximation, E_EXC = 0 and G* = G. The free energy profiles are equivalent to purely thermal reactions.

At Level 0+, E_EXC = E_g, where E_g is the Kohn–Sham band gap of the semiconductor catalyst support. This is a ground-state property accessible through the standard self-consistent field (SCF) iterative procedure and does not involve true excitations of any form. The E_g value can be used as an estimate of the absorption edge and to screen potential semiconductor surfaces for use in photocatalysis. We include it in our analysis to determine the extent of agreement between this and higher-order approximations. For the purposes of this study, we do not consider any bands of Co-d character that appear within the band gap of g-C₃N₄ (midgap states). This omission stems from the fact that in TM-SAC photocatalysts, the principal semiconductor band gap excitation is the major contributor.¹ Metal sites add the possibility of some visible light absorption, but excitations of their localized and short-lived states do not offer useful pathways for efficient electron transfer. This is particularly the case at low loading densities typical of TM-SACs, where the photocatalytic activity is already maximized and further loading (which increases aggregation, metal use, and cost) does not improve the efficiency of the catalyst.¹¹

Level 1, the next step of our protocol, considers E_EXC = VEE, describing an electronic transition that obeys the Franck–Condon rule. In our workflow, charge transfer to the metal site is assumed to be an ultrafast process (≤1 ps for g-C₃N₄ SACs),^12,13 outpacing competing recombination channels. Additionally, it is assumed that the chemical reaction step occurs within this timeframe without geometric relaxation. At this level, the transfer step is described as concerted PCET, in which charge transfer and nuclear rearrangement occur together. The subsequent intermediates relax only to their ground states between the excitations. ΔSCF calculations are performed to obtain VEE values. Without moving from the ground-state optimized geometry, an electron is promoted from the Kohn–Sham VBM of g-C₃N₄ to its CBM. We keep metal midgap states, if any, with their original occupations. The orbitals are allowed to relax self-consistently, with the constrained distribution of electrons across orbitals maintained throughout the procedure. After the calculation with ΔSCF occupations is finished, we compare the calculated total energies to obtain , where E₀ is the total ground-state energy and is the vertical (i.e., without atomic relaxation) excited state energy. VEE is expected to be lower than the band gap E_g, because it corresponds to the optical band gap. The electron and the hole interact through the self-consistent field, reducing the excitation energy relative to the single-particle gap, lowering VEE with respect to the band gap.

At the final step of our protocol, Level 2, , where (rel refers to atomic relaxation) is the excited-state energy with the atomic positions relaxed with non-Aufbau occupations. The underlying physics describes the process of stepwise PCET. The intermediate first relaxes on the photo-excited energy surface before the proton transfer occurs. We expect AEE to always be lower than both E_g and VEE. We note the possibility of considering the excited state thermochemistry (Level 2b) for reactions with a reduced reaction kinetics vs. relaxation of the excitation. This is, however, beyond the scope of the current work in which we consider a non-equilibrium process^14,15 with ground-state thermochemistry (Level 2a).

Among many approaches used to model photochemical processes, ΔSCF stands out as a promising method for screening-level frameworks capable of high-throughput analysis of heterogeneous photocatalysts. This simple approach, the effectiveness of which we demonstrate here, should in the future supplement the existing ab initio molecular dynamics methods, as documented e.g. in recent studies by Carter et al.^16,17 The current approach draws from the existing techniques and goes beyond their typical range of applicability.

To validate our computational protocol, we consider a four-electron OER reaction, going through classical intermediates * + H₂O → *OH → *O → *OOH → *OO → * + O₂ (4e⁻) and a two-electron H₂O₂ evolution * + H₂O → *OH → *(OH)₂ → * + H₂O₂ (2e⁻), which follows the same initial steps but branches after *OH. We allow both 2e⁻ and 4e⁻ oxidation reactions to proceed through the nonclassical *(OH)₂ intermediate.¹⁸ In most of our calculated profiles, O₂ evolution proceeds through the path of *(OH)₂, as the classical *O has a higher ΔG value. We also include the *OO intermediate with the purpose of a more direct comparison of steps that involve the desorption of the final products, *(OH)₂ → * + H₂O₂, and *OO → * + O₂, separating them from the proton-electron transfer steps.

We first turn our attention to the values computed at different approximations for E_EXC by two electron density functionals: a commonly used semilocal PBE+D3(BJ) functional (hereafter PBE) and a hybrid HSE06+D3(BJ) (hereafter HSE), in which vdW D3 dispersion corrections are considered. The PBE results are found to be volatile and without a clear trend that can be easily interpreted (see Fig. S1). Thus, we will focus on the results for the HSE functional, plotted in Fig. 2. In our spin-polarized calculations, we considered two possible excitations for each intermediate, spin-up and spin-down, and selected the one that produced the lowest E_EXC. In the initial calculation without adsorbates, Co has a +3 oxidation state, evidenced by the total magnetic moment of 2.2µ_B for PBE and 2.5µ_B for HSE, depicting the Co³⁺ triplet state, obtained in an unrestricted calculation.

Fig. 2 E_EXC values computed using HSE06+D3(BJ) in different reaction steps: Kohn–Sham band gap (E_g), vertical excitation energy (VEE), and adiabatic excitation energy (AEE).

Out of the three quantities considered as E_EXC, E_g is always the highest, followed by VEE, and finally AEE. The difference between these parameters in subsequent steps is what determines the differences between the calculated ΔG profiles. The VEE to AEE difference for the *OH, *O, *(OH)₂, and *OOH intermediates ranges from low to moderate (0.02–0.10 eV). However, the lowering of AEE with respect to VEE differs substantially for * (bare catalytic site) (1.15 eV) and *OO (0.38 eV). The first leads to a systematic upward shift of the ΔG profile starting from *OH. The second manifests itself with the lowest desorption energy of O₂ found at the HSE-L2a level of just 0.81 eV. On the other hand, the difference between E_g and VEE is less systematic and significant at every step. The biggest changes occur at the second electron transfer steps, *O and *(OH)₂.

Gibbs free energy profiles in the ground state and at higher levels of our protocol are plotted in Fig. S2. We note that the values of absolute Gibbs free energies, especially for gaseous intermediates, depend on the accuracy of the density functional. However, the ΔG profiles considered here rely on relative differences between reaction steps, with the functional-specific errors not influencing the qualitative outcome. Additionally, our calculations of gaseous species are in good agreement with experiments, as shown in Table S3. At Level 0, the 2e⁻ reaction proceeds with *(OH)₂ coupling and desorption as the potential-determining step (PDS) of H₂O₂ production. The 4e⁻ reaction also proceeds through *(OH)₂, as *O is higher in energy by 0.95 eV. In the 4e⁻ process, the coupling to *OOH has the highest PDS energetic cost of 1.97 eV, which is significantly lower than the 2.56 eV for H₂O₂ production. The ground-state calculation thus predicts O₂ to be the majority product over H₂O₂. In the excited state, the potential determining step remains the one involving O–O coupling but depending on the protocol level, its relative intensity changes and influences the predicted major product. Thus, at Level 0+ and 1, *(OH)₂ is more stabilized, and *OOH remains at a similar position in the ΔG diagram. This increases both PDSs and preserves the original preference towards O₂. What changes at Level 2a is a less stable ΔG value of *OH and *(OH)₂, leading to a lower potential difference for HO*–*OH coupling and desorption. At the same time, *OOH is less stable than at lower levels, leading to a lower PDS for the 2e⁻ path towards H₂O₂ over the 4e⁻ path to O₂.

The trends for O₂ vs. H₂O₂ preference are visualized in Fig. 3, including both functionals, and both OER and ORR processes. With the PBE functional, a balance between the formation of the two products is maintained as the points lie close to the X = Y line, indicating equal O₂ and H₂O₂ formation probabilities. Level 2a with PBE (PBE-L2a) is an exception, showing the strongest preference towards H₂O₂, both in 2e⁻ water oxidation (Fig. 3a) and in 2e⁻ O₂ reduction (Fig. 3b) processes. We quantify the magnitude of this preference by introducing the parameter ΔΔG_OER = ΔG_PDS−O2 − ΔG_PDS−H2O2 = 0.98 eV. In the ORR, the thermodynamic feasibility of H₂O₂ production is expressed by G(* + H₂O₂) − G(*OOH); if this value is negative, H₂O₂ is the major product. From Level 0 to Level 1, HSE simulations show a strong preference towards 4e⁻ O₂/H₂O production. At the most precise protocol level, HSE-L2a, both oxidation and reduction paths can be seen as viable to obtain H₂O₂ (ΔΔG = 0.35 eV, G(* + H₂O₂) − G(*OOH)= −0.35 eV), in agreement with existing experimental evidence.^2–4 It is worth noting that the non-classical *(OH)₂ intermediate makes a critical difference in the predicted products. If not present, the OER path could only go through the classical *O intermediate, and H₂O₂ would have only been preferred at the PBE-L2a level (very weakly, ΔΔG = 0.09 eV). However, since the *(OH)₂ reaction path is enabled, H₂O₂ is preferred when PBE-L0 (very weakly, ΔΔG = 0.07 eV), PBE-L2a (ΔΔG = 0.98 eV), or HSE-L2a (ΔΔG = 0.35 eV) are used. While PBE-L0 is qualitatively correct when taken at face value, we attribute the result to fortuitous compensation of systematic errors in the semilocal GGA-type (including PBE) density functionals.¹⁹ This may also be the case in the comparison of HSE-L2a versus PBE-L2a, which gives the same qualitative conclusions in terms of ΔG and ΔΔG values.

Fig. 3 Preference towards different OER and ORR products at different framework levels with the PBE+D3(BJ) and HSE06+D3(BJ) functionals. (a) In the analysis of H₂O oxidation to O₂/H₂O₂, the X and Y axes correspond to the potential determining steps of O₂ and H₂O₂ evolution reactions. X = Y is the line of equal ΔG for potential determining steps. Simulations in which H₂O₂ is preferred are below the line, and O₂ is preferred above the line. (b) For O₂ reduction, Y axis is the Gibbs free energy of the H₂O₂ product versus *OOH intermediate. Below Y = 0, H₂O₂ formation is thermodynamically stable (“downhill”) from *OOH and above, only H₂O is thermodynamically favorable. Orange and purple colors correspond to the 4e⁻ O₂ evolution/4e⁻ O₂ reduction and the 2e⁻ H₂O₂ evolution/2e⁻ O₂ reduction pathways, respectively.

The protocol levels with intermediate accuracy at Levels 0+ and 1 are in agreement in predicting the opposite result. Finally, Level 2a, with both PBE and HSE functionals, gives the correct majority product. The ORR pathways, on the other hand, are not affected by this intermediate, as H₂O₂ production from O₂ does not proceed beyond the 2e⁻ process. The same three functional-protocol level combinations predict H₂O₂ as a thermodynamically allowed ORR product. More specifically, it is determined by the negative ΔG value for the step *OOH + H⁺ + e⁻ → * + H₂O₂ (“downhill” process, purple area of Fig. 3b), in contrast to the positive value that predicts a positive reaction (“uphill”, orange area of Fig. 3b), with only H₂O as a thermodynamically allowed product at lower protocol levels.

In summary, we demonstrated a possible application of a new practical protocol based on computational techniques that are easily transferrable between ground- and excited-state descriptions, while sacrificing little of the efficiency. We showed that, in the case of the OER and ORR on Co/g-C₃N₄, the ΔSCF approach used to relax the excited state (Level 2a) gave predictions of reaction thermodynamics exactly in line with the knowledge from previous experimental studies, with H₂O₂ preference in both reaction paths.

We presented herein a novel, detailed computational protocol to obtain robust free energy profiles of photocatalytic reactions explicitly considering excited state contributions. By building on the ground-state DFT methodology, the ΔSCF treatment of excited-state band occupations allows access to standard methodology and post-processing techniques normally applied to ground-state DFT results.

At this point, we conclude that for the OER and the ORR on Co/g-C₃N₄, the use of a hybrid functional and at least the ΔSCF treatment of the excited states with atomic relaxation is essential to predict the correct products. ΔSCF without relaxation may be insufficient, just as any approximation based on Kohn–Sham eigenvalues. In due course, the protocol should be further validated with calculations for other catalysts and reactions. The models applied here in the minimal demonstration can be expanded to describe interactions with the environment, predict reaction kinetics (in addition to thermodynamics), extend their applicability to other types of reactions (such as photo-electrocatalysis), and, in general, better reproduce experimental conditions.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cc06211j.

Acknowledgements

We gratefully acknowledge Polish high-performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH) for providing computer facilities and support within computational grant no. PLG/2025/017974. We are grateful to the National Science Centre, Poland, for funding (grant no. UMO/2020/39/I/ST4/01446 and UMO-2023/50/E/ST4/00197).

References

1. S. M. Wu and P. Schmuki, Adv. Mater., 2025, 37, e2414889.
2. C. Chu, Q. Zhu, Z. Pan, S. Gupta, D. Huang, Y. Du, S. Weon, Y. Wu, C. Muhich, E. Stavitski, K. Domen and J. H. Kim, Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 6376–6382.
3. W. Wang, Q. Song, Q. Luo, L. Li, X. Huo, S. Chen, J. Li, Y. Li, S. Shi, Y. Yuan, X. Du, K. Zhang and N. Wang, Nat. Commun., 2023, 14, 2493.
4. Y. Shen, Y. Yao, C. Zhu, J. Wu, L. Chen, Q. Fang and S. Song, Chem. Eng. J., 2023, 475, 146383.
5. P. Ren, T. Zhang, N. Jain, H. Y. V. Ching, A. Jaworski, G. Barcaro, S. Monti, J. Silvestre-Albero, V. Celorrio, L. Chouhan, A. Rokicinska, E. Debroye, P. Kustrowski, S. Van Doorslaer, S. Van Aert, S. Bals and S. Das, J. Am. Chem. Soc., 2023, 145, 16584–16596.
6. X. Zhang, H. Su, P. Cui, Y. Cao, Z. Teng, Q. Zhang, Y. Wang, Y. Feng, R. Feng, J. Hou, X. Zhou, P. Ma, H. Hu, K. Wang, C. Wang, L. Gan, Y. Zhao, Q. Liu, T. Zhang and K. Zheng, Nat. Commun., 2023, 14, 7115.
7. C. Feng, J. Luo, C. Chen, S. Zuo, Y. Ren, Z.-P. Wu, M. Hu, S. Ould-Chikh, J. Ruiz-Martínez, Y. Han and H. Zhang, Energy Environ. Sci., 2024, 17, 1520–1530.
8. P. Liu, T. Liang, Y. Li, Z. Zhang, Z. Li, J. Bian and L. Jing, Nat. Commun., 2024, 15, 9224.
9. G. Gao, Y. Jiao, E. R. Waclawik and A. Du, J. Am. Chem. Soc., 2016, 138, 6292–6297.
10. J. M. Herbert, Density-functional theory for electronic excited states, in Theoretical and Computational Photochemistry, ed. G.-I. Cristina and M. Marco, Elsevier, 2023, pp. 69–118.
11. S. Qin, J. Will, H. Kim, N. Denisov, S. Carl, E. Spiecker and P. Schmuki, ACS Energy Lett., 2023, 8, 1209–1214.
12. R. Godin, Y. Wang, M. A. Zwijnenburg, J. Tang and J. R. Durrant, J. Am. Chem. Soc., 2017, 139, 5216–5224.
13. X. Tong, J. Shou, H. Song, Y. Wang, L. Huang and L. Yin, Energy Fuels, 2022, 36, 11532–11541.
14. B. Liu, X. Zhao, C. Terashima, A. Fujishima and K. Nakata, Phys. Chem. Chem. Phys., 2014, 16, 8751–8760.
15. C. Zhang, G. Chen, Y. Si and M. Liu, Phys. Chem. Chem. Phys., 2022, 24, 1237–1261.
16. S. Xu and E. A. Carter, Chem. Rev., 2019, 119, 6631–6669.
17. Y. Yuan, L. Zhou, H. Robatjazi, J. L. Bao, J. Zhou, A. Bayles, L. Yuan, M. Lou, M. Lou, S. Khatiwada, E. A. Carter, P. Nordlander and N. J. Halas, Science, 2022, 378, 889–893.
18. I. Barlocco, L. A. Cipriano, G. Di Liberto and G. Pacchioni, J. Catal., 2023, 417, 351–359.
19. J. P. Perdew, J. Sun, R. M. Martin and B. Delley, Int. J. Quantum Chem., 2016, 116, 847–851.

---