Hydrogen peroxide synthesis on porous graphitic carbon nitride using water as a hydrogen source

Abstract

Using water as a hydrogen source is a promising strategy for alternative hydrogen peroxide (H₂O₂) synthesis. By a series of ab initio molecular dynamics (AIMD) simulations and reactive molecular dynamics (RxMD) calculations, fundamental details have been revealed regarding how liquid water interacts with oxygen on a metal-free carbon nitride catalyst, and the two-step reaction mechanism of H₂O₂ synthesis. Metal-free porous graphitic carbon nitride (g-C₅N₂) catalysts are also systematically screened by using a thermodynamics approach through the ab initio density functional theory (DFT) method. Key results include: (a) pristine g-C₅N₂ is most active to catalyze the H₂O/O₂ reaction and produce H₂O₂; (b) the adsorption and activation of water at unsaturated carbon sites of g-C₅N₂ are critical to initiate the H₂O/O₂ reaction, producing HOO* intermediates; (c) interfacial free water and adsorbed water at g-C₅N₂ form a synergetic proton transfer cluster to promote HOO* intermediates to form H₂O₂. To the best of our knowledge, this work presents long-needed theoretical details of direct H₂O₂ synthesis via the water/oxygen system, which can guide further optimization of carbon-based catalysts for oxygen reduction reactions.

---

1. Introduction

Hydrogen peroxide (H₂O₂) has been extensively used for pulp bleaching, wastewater treatment and green oxidization for chemical synthesis.^1–3 With the increasing global demand, how to synthesize H₂O₂via economic and environment friendly processes is a pressing topic. Currently, over 95% of H₂O₂ is produced using an indirect process, involving energy-intensive multistep anthraquinone oxidation and reduction reactions.^1,4 The industrialized process requires a complex and large-scale infrastructure and causes severe pollution to the environment due to the utilization of aromatic chemicals.⁵ Since the last decade, direct synthesis of H₂O₂via a two-proton hydrogen/oxygen reaction has attracted much attention.^6,7 Various metal and metal oxide catalysts have been proposed and optimized, among which Pd based alloys and nanoparticles seem to be most promising.^8–11 Unfortunately, this method is limited by the heavy use of strong acids and halides, which would promote metal catalyst leaching and as-produced H₂O₂ requires further purification.⁸ In addition, the hydrogen explosion risk at the operational high pressure continues to be a safety hazard and prevents its commercial scale implementation.¹²

Direct H₂O₂ synthesis under ambient conditions, utilizing water as a hydrogen source and combining electrocatalysis or photocatalysis techniques, has witnessed tremendous research efforts recently.^13–17 Kato et al.¹⁸ prepared an Fe–Ru bifunctional catalyst to produce H₂O₂ from H₂O and O₂via visible-light photocatalytic reactions. They observed that H₂O adsorption on Fe sites is critical to the following O₂/H₂O reaction. Electrochemically, H₂O₂ synthesis from H₂O and O₂ is a typical two-proton/two-electron (2H⁺/2e⁻) oxygen reduction reaction (ORR).¹⁹ A number of catalysts have been developed with promising efficiency for electrochemical ORR processes, including noble metal-based electrocatalysts (Au,²⁰ Pt²¹ and Pd¹⁷), single-atom catalysts (Pt@TiN²²), metal oxide catalysts (Fe₃O₄ (ref. 23) and Mn–Ru oxide²⁴) and carbon-based electrocatalysts (N-doped, or B, N co-doped mesoporous carbon^25–28).

Carbon-based catalysts have been considered as an efficient low-cost, metal-free alternative for green and renewable processes.^29–33 For H₂O₂ synthesis from H₂O and O₂, Cui et al.³⁴ reported that oxidized carbon nanotubes (CNTs) exhibited a higher selectivity and a better activity for the two-electron oxygen reduction reaction. Yang and co-workers³⁵ found that epoxy and ring ether groups of graphene oxide exhibit an outstanding electrochemical HO₂⁻ production, achieving a good activity (overpotential <10 mV), an excellent selectivity(≈100%) and a satisfactory stability (over 15 h at 0.45 V in alkaline media). Despite this encouraging progress of carbon-based catalysts for direct H₂O₂ synthesis, detailed fundamental understandings are still incomplete. To just name a few, what is the exact catalytic role of carbon-based materials? what is the critical initial step of H₂O₂ synthesis, water adsorption or oxygen interaction with the catalyst? Is it a one-step reaction or a subsequent two-step proton/electron reaction mechanism? what is the role of the liquid/solid interface?

In this work we report a computational study to reveal the reaction mechanism of direct synthesis of H₂O₂ by a H₂O/O₂ reaction on porous graphitic carbon nitride (g-C₅N₂). Also known as an Aza-fused π-conjugated microporous polymer (Aza-CMP), g-C₅N₂ has a large number of pyridinic nitrogen dopants at zigzag edges, a large surface area, a high pore/volume ratio, and a high electrical conductivity.^36–38 Through a series of AIMD, RxMD and ab initio DFT calculations, the following key results have been elucidated: (a) pristine and hydrogenated g-C₅N₂ catalysts have been screened and the most effective catalyst is partially hydrogenated metastable g-C₅N₂; (b) positively charged carbon sites preferentially chemisorb water molecules, which is critical to promote the H₂O/O₂ interaction and generate HOO* intermediates towards H₂O₂ production; (c) a collection of near-surface water molecules could form a proton transfer chain, thus conveniently promoting the reaction of HOO* intermediates to form H₂O₂.

2. Models and simulation methods

2.1 ReaxFF based reactive molecular dynamics simulation

Developed by van Duin and co-workers, the ReaxFF force field provides affordable and accurate atomic details of reactive complex systems.^39,40 With a bond order concept, ReaxFF force field parameters have been fitted to reproduce ab initio quantum mechanics calculations. Therefore, ReaxFF based reactive molecular dynamics (RxMD) simulations can describe reaction systems with a similar accuracy to ab initio quantum mechanics methods, and effectively handle systems of hundreds of thousands of molecules. The C/N/O/H ReaxFF force field parameters in this work are extracted from the literature.⁴¹ The structural information of g-C₅N₂ has been obtained by both RxMD and ab initio DFT calculations. As presented in Fig S1,† the results of characteristic bonds and angles agree quantitatively well with each other. It is also worth noting that RxMD simulations were performed where only 500 water and 25 oxygen molecules were placed in the simulation box at 300.0 K. The simulation results indicated that this process was only a simple physical mix, and no OOH* and H₂O₂ can be observed after 5 ns (Fig. S2†). RxMD simulations were performed with the LAMMPS software package.⁴² A canonical ensemble was applied where the number of molecules (N), the volume (V), and the temperature (T) remained constant throughout the calculations. The temperature (300.0 K) was maintained by the Nosé–Hoover method with a damping coefficient of 100.0 fs.^43,44 The initial velocities of water and oxygen molecules were assigned according to the Boltzmann distribution. g-C₅N₂ was treated as a flexible solid substrate during the calculation. Newton's equation was integrated by the velocity Verlet algorithm with a time step of 0.25 fs.⁴⁵ A bond-order change of 20% of the original bond length was used to identify the connectivity and molecular species and monitor the system evolution as a function of simulation time.

As illustrated in Fig. 1a, a 2 × 2 supercell was constructed based on the optimized primitive unit cell of g-C₅N₂ containing 30C and 12 N atoms. The simulation box was composed of a single sheet of g-C₅N₂, 68.98 Å (y) × 57.71 Å (x), placed in the middle, and a mixture of 500 water molecules and 25 oxygen molecules, as shown in Fig. 1b. The z dimension of the simulation box was 30.0 Å, so that the interaction between neighboring g-C₅N₂ sheets is negligible. Periodic boundary conditions were applied along the x, y, and z directions. For each RxMD simulation, the system (water, oxygen and g-C₅N₂) was firstly relaxed to optimize the structures, followed by a 2.0 ns calculation to further equilibrate the system. After that, the data were collected for 200 ps for analysis. It is worth noting that there was no restriction on the system: all g-C₅N₂, water and oxygen molecules were allowed to relax and move during the simulation.

Fig. 1 The illustrations of (a) the initial configurations of a 2 × 2 g-C₅N₂ supercell. (b) The simulation box of g-C₅N₂, H₂O and O₂ molecules.

2.2 Ab initio DFT and ab initio MD calculations

Three sets of ab initio quantum mechanics calculations have been performed to: (I) screen and evaluate the structural stability of pristine and hydrogenated g-C₅N₂ models; (II) search and reveal H₂O/O₂ reaction mechanism details, such as the transition state, the reaction pathway and the activation energy barrier; (III) generate the reaction trajectory and dynamics properties, to compare with RxMD calculation results, and to identify the two-proton/two-electron (2H⁺/2e⁻) reaction process.

Geometry optimization and transition state calculation were performed by the ab initio DFT method via the Vienna ab initio simulation package (VASP).^46–48 The exchange-correlation functional was treated via the generalized gradient approximation (GGA) method with Perdew–Burke–Ernzerhof (PBE).⁴⁹ Empirical D3 Grimme's method (DFT-D3) was adopted to improve the description of van der Waals interactions in all calculations.⁵⁰ The planewave kinetic energy cutoff was 450.0 eV. Geometries were optimized until the residual forces were smaller than 0.05 eV Å⁻¹. A 5 × 5 × 1 grid was used for k-point sampling, according to the Monkhorst–Pack method.⁵¹ A vacuum of 20.0 Å was added to the z direction to avoid mirror image interactions. The climbing image-modified nudged elastic band (CI-NEB) method was used to evaluate the activation energies of different reaction paths.⁵² A total of 55 g-C₅N₂ models were generated to represent pristine and all possible hydrogenations in terms of coverage and distribution of hydrogen atoms. Optimized structures and the corresponding equilibrium energies are summarized in Fig. S3 and Table S1 of the ESI.†

Ab initio molecular dynamics (AIMD) calculations were performed by using the CP2K package.^49,53 The system (g-C₅N₂, water and oxygen) was maintained at 300.0 K using the canonical (NVT) ensemble, where the temperature was controlled via the Nosé–Hoover thermostat.^43,54 Each calculation was performed for 2.0 ps with a time step of 0.25 fs. The wave functions were expanded in an optimized double-ζ Gaussian basis set.⁵⁵ The electrostatic energy cutoff for an auxiliary plane wave basis set was 360.0 Ry.⁵⁶ Van der Waals interactions were corrected by the Grimme algorithm (DFT-D3).⁵⁰

2.3 Choice of three simulation methods

While the overall goal of this work is to reveal how water and oxygen react on a carbon nitride catalyst (g-C₅N₂) to produce hydrogen peroxide, the thermal stability of the catalyst and the diffusion/reaction coupled complex process are worth exploring. Thus, we take advantage of three adopted simulation methods to understand the system. First of all, to mimic experimental g-C₅N₂ and provide insight into catalyst screening and optimization, a series of ab initio DFT calculations were performed to calculate the Gibbs free energy and evaluate the thermal stability of all g-C₅N₂ models. Secondly, it is critical to understand the first step of the reaction, the energy barrier and the reaction pathway. Thus, ab initio quantum mechanics methods are appropriate choices. Along this direction, we have performed ab initio DFT and AIMD calculations. Thirdly, H₂O₂ synthesis is a complex process and it occurs at the gas/liquid/solid interface. Obviously, the adsorption and diffusion of reactants (H₂O and O₂) to the solid interface (g-C₅N₂) could be significant parameters. Additionally, after the synthesis reaction, the competitive desorption of H₂O₂ from the solid interface is also vital to a successful separation and purification of the H₂O₂ product. The ReaxFF potential has become a powerful computational tool to explore the chemisorption and reactions of complex systems. More importantly, RxMD simulations not only can describe systems composed of hundreds of thousands of atoms but also can archive a trajectory to the time sale of microseconds. Therefore, RxMD was adopted to study such a complex and multiscale system. The three methods complement each other and provide insights from different perspectives.

3. Results and discussion

3.1 Gibbs free energy calculation and thermal stability of g-C₅N₂ models

The stability of 55 g-C₅N₂ models was firstly investigated by the ab initio DFT based thermodynamics approach.^57–59 The most stable model generally has the lowest Gibbs free energy. In this work, the surface energy (Ω) is calculated according to eqn (1):^60,61where A and G are the surface area and Gibbs free energy of the g-C₅N₂ model, respectively; μ_i is the overall chemical potential of a total number of N_i atoms or molecules of species i. T is the temperature; P is the hydrogen partial pressure. For the solid phase, the change of P has a negligible effect on the Ω value. Thus, eqn (1) could be expressed as:where E_g-C5N2H and E_g-C5N2 are respectively the total energy of hydrogenated g-C₅N₂ and pristine g-C₅N₂ from ab initio DFT calculations; N_H is the number of hydrogen atoms of the g-C₅N₂ model; μ_H is the chemical potential of a single hydrogen atom, which is approximated by the chemical potential of a gas-phase hydrogen molecule:

At the studied temperature of 300.0 K, the gas-phase hydrogen could be treated as an ideal gas. Therefore, the chemical potential could be calculated by using eqn (4):

where is the chemical potential of the hydrogen molecule in the standard state (T⁰ = 298.15 K, P⁰ = 1 atm), is the pressure-dependent term. is the temperature-dependent term, referring to the Gibbs free energy change when the temperature changes from T⁰ = 298.15 K to T = 300.0 K (see Table S2†), while keeping the pressure at P⁰:

The chemical potential of the hydrogen molecule under standard state conditions is:

where is obtained from the NIST thermodynamic tables,⁶² and is calculated using experimental data of heat of formation of water:where the enthalpies of and are adopted from the literature. The oxygen gas entropy in the standard state is from experimental data. Hence, eqn (2) and (3) yield:

As shown in Fig. S3 and Table S1,† the proposed models are constructed based on the degree of hydrogenation, that is, the number of hydrogen atoms added to pristine g-C₅N₂. When multiple candidates exist for a same degree of hydrogenation, the most stable model (with the lowest Gibbs free energy) is then selected as the representative one. With this information, 19 out of the total 55 possible g-C₅N₂ models were selected for the Gibbs free energy calculation to evaluate their thermal stability. The calculations were performed under atmospheric conditions, as a function of variation of hydrogen chemical potential from −0.58 eV to −1.10 eV, which corresponds to the temperature range of 300.0 to 1100.0 K. As shown in Fig. 2a, the Gibbs free energy was expressed as a function of temperature at a fixed partial pressure of hydrogen, P_H2 = 1 atm. A more negative Gibbs free energy indicates a better thermal stability of the g-C₅N₂ model. In addition, the negative slope suggests that the thermal stability generally declines when the temperature increases. According to the Gibbs free energy calculation, the hydrogen coverage (i.e., the degree of hydrogenation) affects the thermal stability. But the dependence is not linear according to the results in Fig. 2a. This is probably due to different activities of C and N sites of g-C₅N₂. The same degree of hydrogenation could have different hydrogen distributions from C and N sites.

Fig. 2 Calculations of the thermodynamic stability and electronic structure. (a) Phase diagrams: the regions of stability of g-C₅N₂ surfaces with different terminations as a function of hydrogen chemical potential variations. (b) Five selected model structures of g-C₅N₂. (c) Electronic localized function (ELF) analysis of the five models. For the scale bar of 0.0 to 1.0, the red region indicates a higher local electron distribution, the green region represents electronic-gas-like pair probability and the blue region shows a higher electronic delocalized distribution.

In this work, we selected five models to study their catalytic roles in H₂O₂ synthesis. As illustrated in Fig. 2b, the models represent respectively pristine (Model 1, no hydrogenation), C-site full hydrogenation (Model 2), N-site full hydrogenation (Model 3), full hydrogenation (Model 4, both C and N sites hydrogenated), and partial hydrogenation (Model 5, C sites fully and N sites partially hydrogenated), which is the most stable model from the Gibbs free energy calculation. The electron localization function (ELF) analysis in Fig. 2c shows the distribution of electrons of the models, where sites with higher electron densities (red color) are preferential to interact with H₂O and O₂ molecules.

3.2 Preferential water adsorption on g-C₅N₂

In order to elucidate which reactant (H₂O or O₂) initiates the reaction, the H₂O/g-C₅N₂ and O₂/g-C₅N₂ systems have been studied separately with the pristine Model 1. Ab initio DFT calculations reveal that O₂ and H₂O adsorption at the pores of g-C₅N₂ has stronger binding energies than that of the basal plane of g-C₅N₂ (Fig. S4†). The adsorption energy is respectively −1.68 eV and −1.14 eV for H₂O/g-C₅N₂ and O₂/g-C₅N₂ at the unsaturated C site, suggesting a preferential interaction between H₂O and g-C₅N₂. The comparison was performed using RxMD simulations where 500H₂O molecules or 50O₂ molecules were placed to interact with pristine g-C₅N₂ (Model 1) at 300.0 K. For the O₂/g-C₅N₂ system, only 4 O₂ molecules were adsorbed at unsaturated C sites. No preferential interaction between O₂/N sites was observed, see Fig. S5 of the ESI.† Meanwhile, for the H₂O/g-C₅N₂ system, as shown in Fig. 3a, water was preferentially distributed to form an ordered ring configuration along the pores of pristine g-C₅N₂. This phenomenon was also observed via ab initio DFT calculation. A zoom-in snapshot of Fig. 3b clearly reveals three types of water molecules in the system, namely, adsorbed water at unsaturated C sites, water molecules near g-C₅N₂ which form a hydrogen bonding network with adsorbed water, and bulk water not shown in Fig. 3b for clarity. The Bader charge analysis in Fig. S6† demonstrates that uncoordinated carbon of pristine g-C₅N₂ carriers a positive charge of 0.18 e by average. This confirms that positively charged carbon of g-C₅N₂ interacts strongly with water through its negatively charged oxygen atom.

Fig. 3 (a) The distribution of water molecules in the pores of dehydrogenated g-C₅N₂. (b) The distribution of water molecules, including adsorbed and bulk H₂O. The atomic density distribution for the water molecule on (c) the C site and (d) the N site. The color codes are: C of g-C₅N₂, gray; N of g-C₅N₂, blue; H of free H₂O, cyan; O of free H₂O, pink; H of adsorbed H₂O, green; O of adsorbed H₂O, yellow. Both line and ball models are used in the snapshots, to emphasize the formation of H₂O₂. For clarity, not all free water molecules are shown.

The distributions of water near C and N sites are shown in Fig. 3c and d, respectively. In Fig. 3c, the first two peaks around 1.50 Å are from absorbed water. As illustrated in Fig. 3b, adsorbed water generally adopts the ‘v’ configuration, atop of C sites. The adsorbed water still demonstrates a certain degree of freedom, which results in two close peaks at 1.50 Å. The peaks around 2.5 Å and 3.7 Å come from hydrogen-bonded water molecules. Limited by the pore size of g-C₅N₂, diameter ∼13.82 Å, only two layers of hydrogen bonds (HB) are allowed, as illustrated by the dotted circles in Fig. 3b. The distribution analysis around N sites revealed one significant peak at around 2.97 Å, which represents hydrogen bonds between two adsorbed water molecules. Due to the negative charge, hydrogen of water stays closer to the N sites, at 2.38 Å.

3.3 Reaction mechanism for direct H₂O₂ synthesis

A series of RxMD and AIMD simulations have been performed to reveal the reaction mechanism of H₂O and O₂ on g-C₅N₂ catalysts and the complete process of direct H₂O₂ synthesis. Fig. 4 presents key snapshots of a representative RxMD simulation, demonstrating how one H₂O₂ molecule was synthesized by a two-step reaction of H₂O/O₂ on pristine g-C₅N₂ (Model 1). It is worth noting that during the RxMD simulation, oxygen and hydrogen were assigned different colors to facilitate tracking of the reaction progress: adsorbed water, oxygen-yellow and hydrogen-green; near surface HB bonded water, oxygen-pink and hydrogen-turquoise; gas-phase oxygen-red. As shown in Fig. 4, step 1 illustrates a critical configuration where water molecules were adsorbed at unsaturated C sites of the pore, and gas-phase oxygen molecules were closer to adsorbed water. From step 1 to step 3, the interaction between O₂ and adsorbed water was witnessed by the OH bond length, changing from 2.037 Å to 1.591 Å, and eventually leading to the dissociation of the adsorbed water and the formation of the OOH* intermediate.

Fig. 4 Proton transfer mechanism (PTM) on the dehydrogenated g-C₅N₂ system: the hydrogen atoms of H₂O₂ come from adsorbed water and bulk water, respectively.

Snapshots from step 4 to step 6 revealed an interesting proton exchange mechanism involving the newly generated OOH* intermediate, a nearby water molecule and another adsorbed water molecule from the neighboring C site. While the overall result is that one oxygen molecule interacts with two adsorbed water molecules to directly synthesize one H₂O₂ molecule, the snapshot of step 6 clearly shows that the two hydrogen atoms of the H₂O₂ molecule came from two different proton donors, which is convincing evidence of the proton exchange event involving near surface water molecules.

It is interesting to note that a different reaction mechanism was also observed from the simulation trajectory, where the O₂ molecule interacted successively with two adsorbed H₂O molecules. As shown by the trajectory snapshots in Fig. 5, the configuration of step 2 is critical: the O₂ molecule diffused to the location where it could simultaneously interact with two adsorbed water molecules. The HOO* intermediate was then produced as a result of the interactions. In addition, since the intermediate was still very close to the other adsorbed water molecule, it could receive the second proton and lead the reaction to completion to produce one H₂O₂ molecule. The color of hydrogen atoms of H₂O₂, step 6 in Fig. 5, clearly demonstrates that no free or hydrogen bonded water participated in the two-step reaction process.

Fig. 5 Direct formation mechanism (DFM) on the dehydrogenated g-C₅N₂ system: the H atoms of H₂O₂ both come from adsorbed H₂O.

The snapshots shown in Fig. 6 are from the AIMD calculation of 5 O₂ molecules, 20H₂O molecules and a pristine g-C₅N₂ model. It is worth noting that, due to the high computational cost, there were fewer oxygen and water molecules in the AIMD calculations. But both RxMD and AIMD were performed at 300.0 K with the same g-C₅N₂ model. The analyses of AIMD and RxMD trajectories give the same conclusion, that is, the first critical step is H₂O adsorption at edge unsaturated C sites, followed by O₂ interaction with adsorbed water to form a HOO* intermediate. Then, different hydrogen suppliers can interact with HOO* to produce H₂O₂. Simultaneously, the interaction between adsorbed water and nearby water molecules was also observed, sequentially producing a hydronium ion (H₃O⁺, step 5) and then a Zundel cation (H₅O₂⁺, step 6). H₂O₂ was eventually synthesized as a reaction product between HOO* and H₅O₂⁺, as illustrated by steps 8 and 9 in Fig. 6. As the reaction proceeded, OH* or O* would be accumulated on the unsaturated edge-C sites. It is probably challenging to regenerate g-C₅N₂ by the H₂O/O₂ system as used in this work. However, convenient techniques are available to recycle the catalyst. Removing OH*/O* and regenerating active sites of carbon-based catalysts can be achieved by electrocatalytic reactions. For example, for 4-electron ORR processes,^63–65 adsorbed OH* is hydrogenated and removed in the form of H₂O. Liang and co-workers studied the 4-electron ORR of nitrogen doped graphene (N-graphene). Their results demonstrated that the removal of adsorbed OH* from the N-graphene surface that is covered by O with a 1/6 monolayer surface oxygen coverage has a small energy barrier of 0.32 eV in the water solution phase.⁶⁶ The regeneration depends closely on the form of deactivated g-C₅N₂ catalysts, the concentration of H₃O⁺ and the applied voltage. To validate the general regeneration process, a proof-of-concept calculation is designed in this work, to mimic a separate electrocatalytic treatment after the successful H₂O₂ production, to remove surface OH and O groups and regenerate the active carbon sites of g-C₅N₂. As shown in Fig. S7(a),† six H₂O and two H₃O⁺ molecules are placed at the pore of the deactivated g-C₅N₂ where its edge nitrogen sites are fully hydrogenated and carbon sites are alternatively bonded with OH and O. The AIMD simulation results in Fig. S7(b) and (c)† show that within 0.2 ps H₃O⁺ interacts with oxygen-containing functional groups at edge carbon sites, transforming OH* and O* back to adsorbed water at edge carbon sites, thus regenerating the g-C₅N₂ catalyst for the next cycle of H₂O₂ production.

Fig. 6 Snapshots of H₂O₂ production on dehydrogenated g-C₅N₂ described by the AIMD simulations. Structures corresponding to the reaction path followed by PTM.

Both RxMD and AIMD reveal the same overall reaction mechanism which can be generally described by two sequential steps: (a) the formation of HOO* intermediates, from the interactions between oxygen and adsorbed water; (b) the completion of H₂O₂ synthesis, by the proton transfer between HOO* and chain-cluster water (or second adsorbed water); (c) H₂O* regeneration, by interactions between residual OH* and H₃O⁺ molecules. The key steps are summarized below:

H₂O adsorption:

H₂O_(bulk) + g-C₅N₂ → g-C₅N₂ *H₂O_(ad) (R1)

HOO* formation:

O_2(gas) + H₂O_(ad) → HOO*_(bulk) + –OH_(dec) (R2)

H₂O₂ formation:

HOO*_(bulk) + H₂O_(chain) → H₂O_2(bulk) + –OH_(dec) (R3.1)

HOO*_(bulk) + *H₂O_(ad) → H₂O_2(bulk) + –OH_(dec) (R3.2)

H₂O chain:

Hydronium ion: *H₂O_(ad) + H₂O_(bulk) → H₃O⁺_(bulk) + –OH_(dec) (R3.1a)

Zundel cation: H₃O⁺_(bulk)+ H₂O_(bulk) → H₅O₂⁺_(bulk) (R3.1b)

H₂O* regeneration: OH* + H₃O⁺_(bulk) → H₂O_(ad) + H₂O_(bulk) (R4)

For other studied g-C₅N₂ models where N- or C-sites are hydrogenated, H₂O₂ was produced by a similar reaction mechanism. The complete process was recorded for a few H₂O₂ molecules, as shown in Fig. S8–11.† Hydrogenated N-sites (Model 3 and Model 4) are potential proton providers. Different H₂O₂ molecules have been identified from the calculations: (a) H₂O₂ synthesized by interacting with an absorbed water and a hydrogenated N-site (Model 3), Fig. S8;† (b) H₂O₂ synthesized by receiving two hydrogens from the same adsorbed water molecule (Model 3), Fig. S9† (c) H₂O₂ synthesized by interacting with an absorbed water and a hydrogenated N-site, via a H₂O proton transfer chain (Model 4), Fig. S10;† (d) H₂O₂ synthesized by directly interacting with two hydrogenated N-sites (Model 4), Fig. S11.† Goclon and Winkler⁶⁷ reported that the oxygen reduction proceeds through a one-step two-electron direct process with a 2.15 eV energy barrier on the amino functionalized g-C₃N₄ structure. In contrast, for the hydrogenated g-C₃N₄ catalyst, the two-step single-electron indirect H₂O₂ synthesis has an energy barrier of 1.78 eV. This indicates that the O₂ molecule can directly interact with hydrogenated sites to produce H₂O₂.

Both RxMD and AIMD calculations confirm the direct H₂O₂ synthesis from chemisorbed water and gaseous colliding oxygen molecules, which is known as the Eley–Rideal (ER) reaction. A direct ER reaction is generally expected to occur only when there is a rather small activation barrier to the reaction, such as a gas-phase radical reactant which undergoes an exothermic reaction. In this work, we adopted the CI-NEB method to study the reaction pathway and activation energy for producing HOO* intermediates and H₂O₂.

As shown in Fig. 7a, for the HOO* intermediate, the initial configuration in the system has a gaseous O₂ and one adsorbed H₂O molecule at g-C₅N₂. While the final configuration has a HOO* intermediate and OH-functionalized g-C₅N₂. For the initial configuration, upon water adsorption at the unsaturated C sites, the O–H bond was elongated, changing from 0.972 Å of the bulk to 1.110 Å. The transition state was identified in which the gaseous O₂ was interacting with the adsorbed water, to the extent that the adsorbed water has two equally elongated O–H bonds, one with the gaseous O₂ and the other from the adsorbed water. The calculation revealed that the activation energy (E_act) for HOO* formation was 0.23 eV, and the overall reaction was exothermic, releasing 0.91 eV from the system. The exothermic nature indicates that the proton transfer from adsorbed water to a nearby O₂ molecule is energetically favored.

Fig. 7 Energy profile and optimized configurations for OOH* and H₂O₂ formation. (a) and (c) Energy profile for OOH* and H₂O₂ formation by DFT calculation, respectively; (b) and (d) the optimized initial state, transition state and final state configurations for OOH* and H₂O₂ formation, respectively.

Fig. 7b shows the reaction pathway by which the HOO* intermediate received the second proton to generate a H₂O₂ molecule. As revealed by RxMD and AIMD calculations, this reaction generally involves ‘free’ water molecules to transfer protons from adsorbed water. It was also observed, see Fig. 5, that the gaseous O₂ could interact successively with two adsorbed H₂O molecules to produce a H₂O₂ molecule, which does not require much diffusion of the HOO* intermediate. It is also worth pointing out that there is no ‘free’ water considered in the CI-NEB calculations, and only the direct reaction between HOO* and adsorbed water to form H₂O₂ was calculated. The result suggested that a larger activation energy (0.34 eV) was necessary to overcome the reaction barrier. But the overall reaction was also exothermic, discharging 0.63 eV from the system.

3.4 Proton transfer along hydrogen bonds

Proton transfer is essential to many aqueous processes, from ion channel function to enzymatic reactions.^68–70 Especially for photocatalytic^71,72 and electrocatalytic reactions,^73,74 simultaneous gain or loss of electrons and protons frequently occur in aqueous solutions. Nevertheless, not enough emphasis has been placed on the effect of proton transfer in aqueous reactions. Experimentally, direct observation of proton transfer is challenging due to the short lifetime of intermediates and the difficulty of distinguishing protons, electrons and their coupling with solvent molecules. For example, Yu and co-workers⁷⁵ investigated the Langmuir–Hinshelwood (LH) and ER mechanisms for the ORR on graphene sheets. They performed ab initio DFT calculations to mimic ER reactions under an alkaline environment, with H shuttling through one or two water molecules. Other theoretical efforts of understanding proton transfer mechanisms have been summarized in recent reviews.^76–85 It has been demonstrated that the proton transfer can largely determine the energy profile and reaction pathway. However, as far as we are aware, there is no quantitative description of proton transfer in ORR processes.

Here, we propose a proton transfer descriptor based on the structural information of involved water molecules. As illustrated in Fig. 8a, the descriptor δ describes the geometric requirement, assuming that proton transfer occurs when the oxygen atoms of water molecules have three identical hydrogen atoms nearby. Therefore, breaking O–H bonds and exchanging hydrogen atoms (protons) require a negligible energy penalty. The δ value is estimated by using the equation:

where θ₀ and L₀ are respectively the equilibrium angle and bond length of free H₂O molecules; θ₁ and θ₂ are the new angles, namely H₁–O–H₃ and H₂–O–H₃ according to Fig. 8a, formed between two nearby water molecules; L₁ is the distance between the oxygen atom and the potential exchangeable H₃ hydrogen (proton). Mathematically, δ reaches a minimum whenever there is a proton transfer event.

Fig. 8 Proton transfer via a chain of water molecules. (a) Schematic representation of the proton transfer process. (b) Trajectory snapshots (I–VII) of critical configurations to the proton transfer for H₂O₂ production. (c) and (d) The change of φ and the total energy in the process of proton transfer, respectively. The dissociated proton exists as a Zundel ion and a hydrated hydronium ion in the systems. C, O, H and N atoms are colored gray, red, white and blue, respectively. The H atom of the adsorbed water that participates in the proton transfer is colored green to guide the viewing.

Previous studies observed Eigen (H₉O₄⁺)^86,87 and H₅O₂⁺ (ref. 88) evolution from bulk water. This suggests that proton transfer events involving two or three water molecules are energetically favorable. In this work, dehydrogenated carbon sites provide localized positive charge centers, which promote H₂O adsorption and subsequent H₅O₂⁺ ion formation (Fig. 8b). Representative snapshots in Fig. 8b demonstrate how near-surface water molecules participate in proton transfer events and promote the HOO* intermediate to receive a proton and thus produce H₂O₂. The AIMD trajectory reveals that two bulk water molecules are involved in this proton transfer process (Fig. 8b(I)): firstly, as shown in Fig. 8b(II), the adsorbed H₂O interacted with a nearby bulk H₂O molecule to form a hydronium H₃O⁺ and an OH functional group of g-C₅N₂, and after that the adsorbed water gave up the proton (green color), see Fig. 8b(III). The hydronium H₃O⁺ then interacted with another bulk water molecule to form a H₅O₂⁺ cation, as shown in Fig. 8b(IV). It is worth noting that the oxygen–oxygen distance, between the second bulk water molecule and the previously formed H₃O⁺, is about 2.70 Å which is within a typical hydrogen bond cutoff.⁸⁹ This suggests that the second water molecule could easily accept a hydrogen atom from H₃O⁺, forming a new hydronium H₃O⁺, see Fig. 8b(V). Since the new H₃O⁺ was much closer to the HOO* intermediate, it easily gives up a proton to HOO* to promote the formation of H₂O₂ (Fig. 8b(VI and VII)). Additionally, direct proton exchange was also observed between adsorbed H₂O molecules and g-C₅N₂, see Fig. S12.†

It has been reported⁸⁹ that a Zundel proton is formed when two water molecules, the oxygen–oxygen distance, are within 1.35 Å. To quantify proton transfer events, we adopted the same cut-off of 1.35 Å to determine whether a proton transfer event could occur between two water molecules. In addition, since the four atoms of the hydronium ion H₃O⁺ adopt a trigonal pyramidal geometry, an angular term is necessary to describe proton transfer events. Using the definition of φ, the dynamics trajectory was monitored. As shown in Fig. 8b, in the process of forming the first H₃O⁺, the value of the bulk and adsorbed H₂O molecules gradually decreases (II to III), eventually reaching the minimum where the bulk water molecule accepted a proton from the adsorbed water, forming H₃O⁺. From (III to IV) of Fig. 8b, the H₃O⁺ interacted with a second water molecule, firstly forming the Zundel proton H₅O₂⁺, where increased to reach the maximum at (IV). The value then decreased and reached the minimum again, where the H₅O₂⁺ decomposed to a new H₃O⁺. The function is similar to the double-minimum potentials that are often used to describe hydrogen bonds and can provide a handy description of proton transfer dynamics. In addition, to better interpret the proton transfer and the H₂O₂ formation, energetics information of each snapshot has been obtained via ab initio DFT. As shown in Fig. 8d, three energy barriers were identified for one successful proton transfer and one HOO hydrogenation event. The energy barrier of forming the first and second H₃O⁺ is 0.20 eV and 0.25 eV (Fig. 8d(I–V)), respectively. More importantly, *HOO only needed to overcome a barrier of 0.24 eV to form H₂O₂ by proton transfer.

3.5 H₂O₂ evolution on different g-C₅N₂ surfaces

The aforementioned results have clearly demonstrated that pristine dehydrogenated g-C₅N₂ is effective to catalyze the H₂O/O₂ reaction and produce H₂O₂. The two-step reaction mechanism depends critically on both adsorbed water (at unsaturated C sites) and free water (near the pores of g-C₅N₂). It is important to evaluate the catalytic performance of other models where C and N sites are partially or fully hydrogenated. Similar to Sec. 2.1, 500H₂O and 25 O₂ molecules were randomly placed on both sides of the g-C₅N₂ model, of the size 68.98 Å (y) × 57.71 Å (x). RxMD calculations were performed at 300.0 K for all five models and the results are shown in Fig. 9a. Firstly, g-C₅N₂ showed a general catalytic activity to H₂O/O₂ reactions and H₂O₂ was produced from all five systems. The best yield (eight H₂O₂) was from the pristine g-C₅N₂ (Model 1, no hydrogenation). And two different types of H₂O₂ molecules were synthesized from different hydrogen resources, see Fig. 9c. This phenomenon was also observed in AIMD simulations of Model 3 (Fig. S8 and S9†) and Model 4 (Fig. S10 and S11†). However, only five H₂O₂ molecules were produced in Model 3 where the C-sites are not hydrogenated and the N sites are fully hydrogenated. The following two observations could explain this result: (a) the transition state calculation indicated that the activation energy of O₂ reaction with H of the edge-N-H species is 0.58 eV (Fig S14(a)†), which is higher than that of O₂ reaction with adsorbed water at the edge C sites. The activation energy for HOO* and H₂O₂ formation was respectively only 0.23 eV and 0.32 eV for Model 1. Thus, O₂ shall preferentially interact with the adsorbed H₂O. (b) As shown in Fig S8,† AIMD calculations reveal that when adsorbed H₂O of edge C sites and hydrogenated H of edge N sites coexist, O₂ preferentially reacted with adsorbed H₂O to form HOO* intermediates. Those HOO* intermediates could then capture H from those edge N sites to form H₂O₂. Based on these observations, we conclude that O₂ reaction with two adsorbed water molecules, as shown for Model 1, is energetically favored and that Model 1 is more active than Mode 3 for the studied H₂O₂ synthesis.

Fig. 9 RxMD simulations of H₂O₂ synthesis on five g-C₅N₂ models at 300.0 K. (a) The number of H₂O₂ molecules produced as a function of simulation time. (b) Snapshots of H₂O and O₂ at the dehydrogenated g-C₅N₂ phase surface described by the RxMD simulations for 200.0 ps. (c) Observed H₂O₂ production from different nanopores of Model 1 as highlighted in (b).

For Model 2 (C-site fully hydrogenated), Model 4 (C-site and N-site fully hydrogenated) and Model 5 (the most stable model, with C sites fully and N sites partially hydrogenated), the yield is relatively lower, only 1, 3 and 2H₂O₂ molecules were produced respectively during the simulation of 1000 ps. According to AIMD results, O₂ can only capture protons from the hydrogenated N sites due to the very stable C–H bonds. In addition, the adsorption configurations of H₂O and O₂ on Model 5 have been optimized via ab initio DFT calculations based on Bader charge analysis (Fig. S13a†). The adsorption energy of O₂ is −0.95 eV (Fig. S13b†), and that of H₂O is only −0.05 eV (Fig. S13c†). Stronger O₂ adsorption can significantly weaken the bond energy of O₂, which will reduce the selectivity of H₂O₂ and produce water as a side reaction product. Thus, the number of formed H₂O₂ molecules of Model 4 are more than that of Model 2 and Model 5 with all the edge-C sites hydrogenated.

Although the Gibbs free energy was calculated to evaluate possible g-C₅N₂ hydrogenation, other factors, such as defects and solutions, could also affect the stability and reactivity of g-C₅N₂ catalysts. If one only considers the best reactivity, as shown in Fig. 9a, the pristine g-C₅N₂, in which neither edge nitrogen nor carbon sites are hydrogenated, demonstrates the best catalytic performance, while the most stable g-C₅N₂ model (Model 5) has a relatively poor catalytic performance towards H₂O₂ synthesis. The activity is the major contributor to H₂O₂ synthesis, but the successful application of g-C₅N₂ materials shall also depend strongly on the stability. Regardless of the mechanism or yield of H₂O₂ formation analysis, both suggest that unsaturated edge carbon sites do not promote immediate water dissociation and could be reactivated afterwards. We expect that those unsaturated edge sites are available both from the synthesized g-C₅N₂ and during the process of catalyzing H₂O/O₂ reactions.

4. Conclusion

In summary, this work represents a computational study of direct H₂O₂ synthesis from H₂O and O₂ on a metal-free porous graphitic carbon nitride catalyst. Through a combination of ab initio molecular dynamics simulation, reactive molecular dynamics (RxMD) calculation, and the ab inito density functional theory based thermodynamics approach, g-C₅N₂ models with different hydrogenation have been examined to catalyze H₂O/O₂ reactions for H₂O₂ production. The Gibbs free energy calculation and the thermal stability discussion of various g-C₅N₂ models could be a general protocol to screen and evaluate carbon-based materials which are doped with a wide range of other elements. The two-step reaction mechanism has been elucidated, which involves sequentially the adsorption and activation of water at unsaturated C sites of g-C₅N₂, producing HOO* intermediates, and synergetic proton transfer to promote HOO*-to-H₂O₂ reactions, via clusters of free and adsorbed water. Other fundamental facts include that pristine g-C₅N₂ is the most active catalyst and that a preferential water adsorption (instead of oxygen adsorption) on g-C₅N₂ is critical to direct H₂O₂ synthesis. Future development of new catalysts could take this as a general criterion to design reactive sites and optimize catalyst geometry, to simultaneously favor water activation at reactive sites and water-chain assisted fast proton transfer processes. This work paves the way for using carbon-based sustainable catalysts for direct H₂O₂ synthesis via this promising water/oxygen strategy. In addition, the fundamental insight could also shed light on selecting and designing other green catalysts for heterogeneous oxygen hydrogenation reactions.

Conflicts of interest

The authors declare no competing financial interests.

Acknowledgements

Y. Y. Cao and J. G. Wang acknowledge the support by the National Natural Science Foundation of China (Grant No 21625604, 21671172, 21776251 and 21706229). L.L. Huang acknowledges the U.S. National Science Foundation (NSF) for support through Grant CHE-1710102. We are very pleased to thank the OU Supercomputing Center for Education & Research (OSCER) at the University of Oklahoma for their dedicated support.

References

1. J. M. Campos-Martin, G. Blanco-Brieva and J. L. G. Fierro, Angew. Chem., Int. Ed., 2006, 45, 6962–6984.
2. E. Brillas, I. Sires and M. A. Oturan, Chem. Rev., 2009, 109, 6570–6631.
3. N. Agarwal, S. J. Freakley, R. U. McVicker, S. M. Althahban and N. Dimitratos, Science, 2017, 358, 223–227.
4. Y. Yi, L. Wang, G. Li and H. Guo, Catal. Sci. Technol., 2016, 6, 1593–1610.
5. J. Garcia-Serna, T. Moreno, P. Biasi, M. J. Cocero, J. P. Mikkola and T. O. Salmi, Green Chem., 2014, 16, 2320–2343.
6. J. K. Edwards, B. Solsona, E. N. Ntainjua, A. F. Carley, A. A. Herzing, C. J. Kiely and G. J. Hutchings, Science, 2009, 323, 1037–1041.
7. Q. Liu, J. C. Bauer, R. E. Schaak and J. H. Lunsford, Angew. Chem., Int. Ed., 2008, 47, 6221–6224.
8. G. M. Lari, B. Puertolas, M. Shahrokhi, N. Lopez and J. Perez-Ramirez, Angew. Chem., Int. Ed., 2017, 56, 1775–1779.
9. J. S. Jirkovský, I. Panas, E. Ahlberg, M. Halasa, S. Romani and D. J. Schiffrin, J. Am. Chem. Soc., 2011, 133, 19432–19441.
10. I. Kim, M.-g. Seo, C. Choi, J. S. Kim, E. Jung, G.-H. Han, J.-C. Lee, S. S. Han, J.-P. Ahn, Y. Jung, K.-Y. Lee and T. Yu, ACS Appl. Mater. Interfaces, 2018, 10, 38109–38116.
11. E. Pizzutilo, S. J. Freakley, S. Cherevko, S. Venkatesan, G. J. Hutchings, C. H. Liebscher, G. Dehm and K. J. J. Mayrhofer, ACS Catal., 2017, 7, 5699–5705.
12. Y. Shiraishi, Y. Kofuji, H. Sakamoto, S. Tanaka, S. Ichikawa and T. Hirai, ACS Catal., 2015, 5, 3058–3066.
13. Y. Isaka, Y. Kawase, Y. Kuwahara, K. Mori and H. Yamashita, Angew. Chem., Int. Ed., 2019, 58, 5402–5406.
14. S. Zhao, T. Guo, X. Li, T. Xu, B. Yang and X. Zhao, Appl. Catal., B, 2018, 224, 725–732.
15. M. Gryszel, A. Markov, M. Vagin and E. D. Głowacki, J. Mater. Chem. A, 2018, 6, 24709–24716.
16. X. Shi, S. Siahrostami, G.-L. Li, Y. Zhang, P. Chakthranont, F. Studt, T. F. Jaramillo, X. Zheng and J. K. Nørskov, Nat. Commun., 2017, 8, 701.
17. S. Siahrostami, A. Verdaguer-Casadevall, M. Karamad, D. Deiana, P. Malacrida, B. Wickman, M. Escudero-Escribano, E. A. Paoli, R. Frydendal, T. W. Hansen, I. Chorkendorff, I. E. L. Stephens and J. Rossmeisl, Nat. Mater., 2013, 12, 1137–1143.
18. S. Kato, J. Jung, T. Suenobu and S. Fukuzumi, Energy Environ. Sci., 2013, 6, 3756.
19. I. Yamanaka and T. Murayama, Angew. Chem., Int. Ed., 2008, 47, 1900–1902.
20. H. Jing, Q. F. Zhang, N. Large, C. M. Yu, D. A. Blom, P. Nordlander and H. Wang, Nano Lett., 2014, 14, 3674–3682.
21. C. H. Choi, H. C. Kwon, S. Yook, H. Shin, H. Kim and M. Choi, J. Phys. Chem. C, 2014, 118, 30063–30070.
22. S. Yang, J. Kim, Y. J. Tak, A. Soon and H. Lee, Angew. Chem., Int. Ed., 2016, 55, 2058–2062.
23. W. R. P. Barros, Q. L. Wei, G. X. Zhang, S. H. Sun, M. R. V. Lanza and A. C. Tavares, Electrochim. Acta, 2015, 162, 263–270.
24. M. B. Zakaria, C. Li, M. Pramanik, Y. Tsujimoto, M. Hu, V. Malgras, S. Tominaka and Y. Yamauchi, J. Mater. Chem. A, 2016, 4, 9266–9274.
25. 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, J. Am. Chem. Soc., 2018, 140, 7851–7859.
26. T.-P. Fellinger, F. Hasché, P. Strasser and M. Antonietti, J. Am. Chem. Soc., 2012, 134, 4072–4075.
27. L. Han, Y. Sun, S. Li, C. Cheng, C. E. Halbig, P. Feicht, J. L. Hübner, P. Strasser and S. Eigler, ACS Catal., 2019, 9, 1283–1288.
28. J. Park, Y. Nabae, T. Hayakawa and M. a. Kakimoto, ACS Catal., 2014, 4, 3749–3754.
29. Y. Cao, S. Deng, Q. Fang, X. Sun, C. Zhao, J. Zheng, Y. Gao, H. Zhuo, Y. Li, Z. Yao, Z. Wei, X. Zhong, G. Zhuang and J. Wang, Nanotechnology, 2019, 30, 335403.
30. K. Gong, F. Du, Z. Xia, M. Durstock and L. Dai, Science, 2009, 323, 760–764.
31. Z. Pan, K. Wang, Y. Wang, P. Tsiakaras and S. Song, Appl. Catal., B, 2018, 237, 392–400.
32. L. Tao, Q. Wang, S. Dou, Z. Ma, J. Huo, S. Wang and L. Dai, Chem. Commun., 2016, 52, 2764–2767.
33. P. Zhang, D. Sun, A. Cho, S. Weon, S. Lee, J. Lee, J. W. Han, D.-P. Kim and W. Choi, Nat. Commun., 2019, 10, 940.
34. Z. Lu, G. Chen, S. Siahrostami, Z. Chen, K. Liu, J. Xie, L. Liao, T. Wu, D. Lin, Y. Liu, T. F. Jaramillo, J. K. Norskov and Y. Cui, Nat. Catal., 2018, 1, 156–162.
35. H. W. Kim, M. B. Ross, N. Kornienko, L. Zhang, J. Guo, P. Yang and B. D. McCloskey, Nat. Catal., 2018, 1, 282–290.
36. V. Briega-Martos, A. Ferre-Vilaplana, A. de la Peña, J. L. Segura, F. Zamora, J. M. Feliu and E. Herrero, ACS Catal., 2016, 7, 1015–1024.
37. Y. Kou, Y. Xu, Z. Guo and D. Jiang, Angew. Chem., Int. Ed., 2011, 50, 8753–8757.
38. Z.-D. Yang, W. Wu and X. C. Zeng, J. Mater. Chem. C, 2014, 2, 2902–2907.
39. A. C. T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard, J. Phys. Chem. A, 2001, 105, 9396–9409.
40. K. Chenoweth, A. C. T. van Duin and W. A. Goddard III, J. Phys. Chem. A, 2008, 112, 1040–1053.
41. L. L. Huang, M. Seredych, T. J. Bandosz, A. C. T. van Duin, X. H. Lu and K. E. Gubbins, J. Chem. Phys., 2013, 139, 194707.
42. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
43. S. Nose, J. Chem. Phys., 1984, 81, 511–519.
44. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
45. L. Verlet, Phys. Rev., 1967, 159, 98–103.
46. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
47. G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50.
48. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
49. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
50. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
51. R. A. Evarestov and V. P. Smirnov, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 233101.
52. G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000, 113, 9901–9904.
53. J. VandeVondele, M. Krack, F. Mohamed, M. Parrinello, T. Chassaing and J. Hutter, Comput. Phys. Commun., 2005, 167, 103–128.
54. W. G. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
55. J. VandeVondele and J. Hutter, J. Chem. Phys., 2007, 127, 114105.
56. G. Lippert, J. Hutter and M. Parrinello, Mol. Phys., 1997, 92, 477–487.
57. F. Bottin, F. Finocchi and C. Noguera, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 13.
58. E. Heifets, S. Piskunov, E. A. Kotomin, Y. F. Zhukovskii and D. E. Ellis, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 13.
59. Q. Cai, J.-g. Wang, Y. Wang and D. Mei, J. Phys. Chem. C, 2016, 120, 19087–19096.
60. D. Gao, H. Zhou, F. Cai, D. Wang, Y. Hu, B. Jiang, W.-B. Cai, X. Chen, R. Si, F. Yang, S. Miao, J. Wang, G. Wang and X. Bao, Nano Res., 2017, 10, 2181–2191.
61. K. Reuter and M. Scheffler, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 68, 045407.
62. http://https://webbook.nist.gov/ .
63. R. Cao, J.-S. Lee, M. Liu and J. Cho, Adv. Energy Mater., 2012, 2, 816–829.
64. Z. L. Wang, D. Xu, J. J. Xu and X. B. Zhang, Chem. Soc. Rev., 2014, 43, 7746–7786.
65. R. Ma, G. Lin, Y. Zhou, Q. Liu, T. Zhang, G. Shan, M. Yang and J. Wang, npj Comput. Mater., 2019, 5, 78.
66. L. Yu, X. Pan, X. Cao, P. Hu and X. Bao, J. Catal., 2011, 282, 183–190.
67. J. Goclon and K. Winkler, Appl. Surf. Sci., 2018, 462, 134–141.
68. O. Gerlits, T. Wymore, A. Das, C. H. Shen, J. M. Parks, J. C. Smith, K. L. Weiss, D. A. Keen, M. P. Blakeley, J. M. Louis, P. Langan, I. T. Weber and A. Kovalevsky, Angew. Chem., Int. Ed., 2016, 55, 4924–4927.
69. C. T. Supuran, Biochem. J., 2016, 473, 2023–2032.
70. E. Nango, A. Royant, M. Kubo, T. Nakane, C. Wickstrand, T. Kimura, T. Tanaka, K. Tono, C. Y. Song, R. Tanaka, T. Arima, A. Yamashita, J. Kobayashi, T. Hosaka, E. Mizohata, P. Nogly, M. Sugahara, D. Nam, T. Nomura, T. Shimamura, D. Im, T. Fujiwara, Y. Yamanaka, B. Jeon, T. Nishizawa, K. Oda, M. Fukuda, R. Andersson, P. Bath, R. Dods, J. Davidsson, S. Matsuoka, S. Kawatake, M. Murata, O. Nureki, S. Owada, T. Kameshima, T. Hatsui, Y. Joti, G. Schertler, M. Yabashi, A. N. Bondar, J. Standfuss, R. Neutze and S. Iwata, Science, 2016, 354, 1552–1557.
71. C. Chen, T. Shi, W. Chang and J. Zhao, ChemCatChem, 2015, 7, 724–731.
72. K. Yamamoto and K. Takatsuka, Phys. Chem. Chem. Phys., 2018, 20, 6708–6725.
73. A. Badalyan and S. S. Stahl, Nature, 2016, 535, 406–410.
74. C. W. Lee, J. S. Hong, K. D. Yang, K. Jin, J. H. Lee, H.-Y. Ahn, H. Seo, N.-E. Sung and K. T. Nam, ACS Catal., 2018, 8, 931–937.
75. Q. Ly, B. V. Merinov, H. Xiao, W. A. Goddard and T. H. Yu, J. Phys. Chem. C, 2017, 121, 24408–24417.
76. M. T. M. Koper, Chem. Sci., 2013, 4, 2710–2723.
77. U. W. Schmitt and G. A. Voth, J. Chem. Phys., 1999, 111, 9361–9381.
78. S. Hammes-Schiffer, Acc. Chem. Res., 2001, 34, 273–281.
79. J. M. Saveant, Acc. Chem. Res., 1993, 26, 455–461.
80. S. Scheiner, Acc. Chem. Res., 1985, 18, 174–180.
81. M. H. V. Huynh and T. J. Meyer, Chem. Rev., 2007, 107, 5004–5064.
82. D. R. Weinberg, C. J. Gagliardi, J. F. Hull, C. F. Murphy, C. A. Kent, B. C. Westlake, A. Paul, D. H. Ess, D. G. McCafferty and T. J. Meyer, Chem. Rev., 2012, 112, 4016–4093.
83. R. I. Cukier and D. G. Nocera, Annu. Rev. Phys. Chem., 1998, 49, 337–369.
84. J. M. Mayer, Annu. Rev. Phys. Chem., 2004, 55, 363–390.
85. S. Hammes-Schiffer and A. A. Stuchebrukhov, Chem. Rev., 2010, 110, 6939–6960.
86. M. Eigen, Angew. Chem., 1963, 75, 489–508.
87. M. Eigen, Angew. Chem., Int. Ed., 1964, 3, 1–72.
88. G. Zundel, Adv. Chem. Phys., 2000, 111, 1–217.
89. D. Marx, ChemPhysChem, 2006, 7, 1848–1870.

---

Footnote

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

---