Mode-selective H₂O dissociation on Pt(111) under two-dimensional confinement

Abstract

Understanding how spatial confinement alters chemical reactivity at surfaces remains an unanswered but fundamental question. In this work, we investigate the dissociative chemisorption of H₂O on Pt(111), focusing on how two-dimensional (2D) confinement, introduced by graphene (Gr), boron nitride (BN), and graphitic carbon nitride (C₃N₄), modulates the mode-selectivity of this reaction. Building on our earlier findings of barrier reduction under such confinements, we employ a reaction path Hamiltonian approach to uncover how vibrational energy redistribution is affected in both vibrationally adiabatic and non-adiabatic regimes. Our findings reveal that 2D confinement significantly enhances dissociation probabilities for both ground state and vibrationally excited state. Remarkably, the ground-state reactivity under confinement reaches levels comparable to those achieved via single-quantum vibrational excitation of H₂O on the bare surface. The symmetric stretching and bending modes have more influence on reactivity owing to their greater softening near the transition state. Non-adiabatic factors mediated by curvature- and Coriolis-couplings influence energy redistribution by amplifying low-energy reactivity. In the unconfined case, Coriolis interactions dominate and facilitate over-the-barrier processes through vibrational de-excitation. Among the confined cases, Gr cover exhibits the highest effect on reactivity, yielding nearly barrierless dissociation pathways through strong coupling between the reaction path and vibrationally excited symmetric and bending modes. Altogether, the analysis of spatial confinement and vibrational excitation in this work paves the way for rational catalytic design.

---

1 Introduction

The quest for high-performance catalysts has long focused on enhancing activity and selectivity, both of which are profoundly influenced by the local environment surrounding the active site.¹ Considerable efforts have been devoted toward designing catalytic systems that strike an optimal balance between reactivity and cost, such as, bimetallic alloy surfaces composed of transition metals (e.g., Cu/Ni, Ni/Pt) with tunable composition,^2–9 that mitigate the limitations of monometallic catalysts. A recent strategy for modulating catalytic efficiency involves leveraging spatial confinement to tailor the chemical environment around active sites. Confined reaction spaces, such as those found in metal–organic frameworks, zeolites, or nanotubes, have enabled altered reaction pathways and stabilization of key intermediates.^10–17 Among the various forms of confinement, planar confinement introduced through the use of atomically thin two-dimensional (2D) materials supported on metal surfaces offers a more tractable platform for probing chemical reactivity under spatial restriction.^18–24 These heterostructures provide well-defined nano-interface that can host reactions under spatial constraints, thereby allowing detailed investigations of catalytic processes in a controlled manner.^25–30 The tunable nature of such systems, permitted through the choice of 2D overlayer, substrate, and stacking configuration, had been harnessed to manipulate the microenvironment around catalytic active sites.^31–33 The ability of 2D overlayers to influence adsorption energetics, manipulate surface charge distributions, and even redirect reaction paths has positioned 2D confinement as a promising approach for enhancing catalytic performance across a range of important chemical transformations, including those relevant in heterogeneous catalysis and electrocatalysis.^34–44

Water (H₂O) dissociation on metal surfaces, especially on Pt-based catalysts, has been extensively investigated owing to its fundamental and industrial significance in processes such as water–gas shift (WGS) reaction.^45–50 Both experimental techniques and density functional theory (DFT) calculations have uncovered the detailed energetics and geometries of H₂O molecules on Pt surface, revealing key insights into coverage-dependent behavior and the role of surface morphology, such as steps and kinks, in lowering activation barriers for dissociation.^51–54 A wealth of studies has established Pt(111) as a prototypical substrate for exploring the mechanisms of H₂O adsorption and dissociation, owing to its ordered crystal structure and high catalytic activity.^51,54,55 Under suitable conditions, H₂O can intercalate beneath graphene (Gr) overlayers on metal substrates,^56,57 introducing an additional level of confinement that can markedly influence surface chemistry. Reactions involving H₂O under such confined environments have therefore attracted growing interest, with reports demonstrating that nanoconfinement can modify the kinetics and thermodynamics of H₂O formation, dissociation, and oxidation.^16,58–60 For instance, experimental studies by Wang et al. revealed that the confined space at the bilayer-silica/Ru(0001) interface accelerates H₂O production,⁶¹ while Xie et al. used porous Gr layers to enhance photoelectrochemical H₂O oxidation performance.⁶² Additionally, investigations on Gr/Ir(111) systems have shown that the encapsulated environment can enhance the stability of the reactive intermediates and influence the hydrogenation of oxygen.⁶³ Complementary to these observations, Politano and co-workers demonstrated, through combined experimental and theoretical analyses, that H₂O can migrate to the Gr/Cu(111) interface and dissociate into OH and H species with a low energy barrier (0.62 eV), accompanied by partial decoupling of Gr edges from the substrate.⁶⁴ While neglecting possible chemical modifications of the 2D overlayers arising from their interaction with H₂O, our prior theoretical investigation using DFT based calculations suggested that confinement by 2D Gr, boron nitride (BN), and graphitic carbon nitride (C₃N₄) significantly impacts the energetic landscape of H₂O dissociation on Pt(111), leading to differing degrees of adsorption weakening and barrier reduction.⁶⁵ These effects arise due to the combination of mechanical, geometric, and electronic influences introduced by the confined environment. Furthermore, by applying the harmonic transition state theory, we evaluated the tunneling-corrected rate constants for the dissociative chemisorption of H₂O,⁶⁵ providing an approximate measure of reactivity under confinement. These findings motivated a deeper investigation into the dynamical intricacies of this reaction. The course of a reaction can be controlled by exciting specific vibrational modes, a phenomenon known as mode selectivity. While confinement effect on activation energies and site accessibility are well-discussed, much less is known about how 2D confinement tunes mode selectivity. Infrared radiation applied to a molecule undergoing dissociation on a metal surface often leads to creation of electron–hole pairs and subsequent electronic excitations. However, several experiments have successfully irradiated the vibrational modes of adsorbate molecules, thereby opening the avenue to realize quantum state resolved reactivity during dissociative chemisorption processes on a metal surfaces. In the past decade, experimentalists have been able to selectively prepare surface-incident H₂O and CH₄ with one or two quanta excitations of a particular vibrational mode via state-resolved infrared (IR) laser pumping in a molecular beam.^66–68 More recently, femtosecond time-resolved surface-enhanced Raman scattering experiments has also enabled the vibrational excitations of adsorbates by utilizing ultrafast laser pulses.⁶⁹ These findings motivated us to investigate the effect of selective mode excitation under the influence of 2D confinement, which could provide theoretical benchmark for future experimental measurements on dissociative chemisorption.

For an accurate analysis of vibrational energy redistribution in mode-specific and bond-selective chemistry, full-dimensional quantum simulations are essential. However, these simulations are computationally demanding. Although full-dimensional potential energy surfaces (PESs) provide chemically accurate descriptions of molecule–surface interaction energies and reaction barriers, their construction and subsequent quantum dynamical treatments pose major challenges. Nonetheless, researchers have made impressive strides by developing multi-dimensional PESs with a large number of degrees-of-freedom (DOFs)^70–74 enabling quantum dynamics studies with increasingly realistic representations and enhanced precision. In the context of H₂O dissociation, Jiang and Guo⁷⁵ constructed a 9D PES for the Ni(111) surface using a permutation-invariant polynomial neural network (PIP-NN),⁷⁶ which facilitated subsequent quantum and classical simulations to examine the role of vibrational and rotational excitation, as well as site-specific adsorption.⁷⁷ These simulations offered a semi-quantitative explanation of mode-specific behavior observed in quantum-state-resolved molecular beam experiments.⁶⁸ Following this, Zhang and coworkers performed 7D^78,79 as well as 9D⁸⁰ quantum dynamical simulations on a high-quality 9D PES, constructed using neural network (NN) approach, to evaluate mode-specific reactivity for H₂O dissociation on Cu(111) surfaces. Furthermore, to investigate the effect of vibrational mode excitation during H₂O dissociation on bare Ni(100) surface, 9D quantum dynamics calculations were carried out by Liu et al.^81,82 on “accurate” 9D potential energy surfaces. Despite these advancements, the computational cost of multi-dimensional quantum dynamical calculations for polyatomic molecules remains prohibitive, necessitating the exploration of alternative approaches.^83–85 One such approach is an approximate dynamical method developed by Jackson and co-workers,^86–91 which is based on the quantum version of the reaction path Hamiltonian (RPH) theory originally proposed by Miller and coworkers.⁹² The reaction path is defined as an abstract 1D path connecting the saddle point to the minima, comprised of a combination of 3N − 6 DOFs. It is assumed that the PES remains harmonic with respect to displacements away from the reaction path. The main computational advantage of this method lies in the need to compute only a one-dimensional path, requiring fewer ab initio points compared to a multidimensional PES. When implemented with the inclusion of various normal modes, this method has proven effective in capturing the correct trends of tunneling probabilities, mode selectivity and surface temperature effects via the curvature effects and surface lattice motions.^86,87,93

The influence of 2D confinement on mode-selective reactivity is largely an uncharted frontier. To address this gap, we apply the RPH-based quantum dynamics for investigating and comparing the dissociative chemisorption of H₂O on Pt(111) surface, with and without confinement by 2D layers: Gr, BN, and C₃N₄. To the best of our knowledge, this is the first effort to elucidate the effect of 2D confinement on vibrationally mediated H₂O dissociation on metal surfaces. Our goal, in this present article, is to scrutinize the potential changes in vibrational energy transfer, mode selectivity, and overall dissociation probability under different confining covers. The dissociation probabilities are evaluated within the vibrationally adiabatic approximation and further refined through a non-adiabatic treatment that captures the couplings among the normal modes and the reaction path.

2 Methodology & computational details

The energetics calculations have been performed using the Vienna ab initio simulation package (Vasp 5.4.1) based on DFT.^94–97 The plane wave basis set expansion is truncated at a kinetic energy cut-off of 400 eV. To treat the exchange–correlation effects, the Perdew–Burke–Erzenhof (PBE) functional with generalized gradient approximation (GGA)^98,99 is used. The DFT-D3 correction with Becke–Johnson damping function^100,101 has been employed to incorporate the dispersion forces between the confining surfaces and H₂O.^102,103 For ionic core–electron interactions, the projector-augmented wave (PAW)^104,105 pseudopotential is used. Two-dimensional confinements have been modeled by using commensurate supercells of Pt and 2D overlayers. Specifically, 3 × 3 graphene, 3 × 3 BN and 1 × 1 C₃N₄ monolayers are placed above a four-layer slab of Pt(111), resulting in lattice mismatches of 0.84%, 1.98%, and 2.83%, respectively, as detailed in our previous work.⁶⁵ These values are sufficiently small to be considered negligible in terms of their effect on the overall reactivity.¹⁰⁶ The periodic images along the z-direction are separated by an optimal vacuum slab to avoid unphysical interferences. The geometries of the initial and final states are optimized with the bottom two layers of the Pt(111) surface fixed. Fig. S1 presents the optimized initial state on Pt(111) as well as the confining environments of the H₂O molecule. The most stable adsorption geometry of an isolated H₂O molecule on Pt(111) surface corresponds to a nearly flat configuration, where the oxygen atom binds through its lone pair to a top Pt site and the molecular plane is almost parallel to the surface.^{51,54,102,107} We adopted this orientation as the initial structure for evaluating the reaction energy profile on bare Pt(111). Spin-polarization is included in all calculations and the reciprocal space is sampled using a 4 × 4 × 1 Gamma (Γ) centered k-point grid. A convergence criterion of 0.01 eV Å⁻¹ is set for the ionic forces.

2.1 Reaction path Hamiltonian approach

Molecular geometries for H₂O dissociation along the reaction path have been determined using the climbing image-nudged elastic band (CI-NEB) method,¹⁰⁸ while keeping the surface atoms fixed. The CI-NEB calculations are considered as converged when the corresponding ionic forces reach below 0.05 eV Å⁻¹.

To construct a reaction path Hamiltonian, one first needs to locate the reaction path, or minimum energy path (MEP), defined usually in terms of mass-weighted cartesian coordinates as the path of steepest descent from a saddle point on the potential energy surface to the respective minima in the reactant and product configurations. After the identification of the reactant (H₂O at the top site on the surface) and product (chemisorbed H and OH) states, subsequent calculations are performed with the CI-NEB method to obtain the TS and the molecular configurations at several points along the reaction path (RP). The distance between successive configurations along the MEP, i.e., the arc length along the reaction path, is represented by s, and is defined as: , where x_i are the mass-weighted Cartesian coordinates of the atoms of the H₂O molecule. The upper limit in this summation represents the total number of DOFs of the concerned molecule, i.e., 3N for an N-atomic molecule. The configurational space in the RPH method is described by the arc length-coordinate (s) along the MEP and the (3N − 7) normal mode vectors orthogonal to this path, obtained from diagonalization of the force constant matrix. The remaining six eigenvectors correspond to negligible translational and rotational motions. The coordinate set representing a reacting system in the (3N − 6)-dimensional configuration space, decomposed into a 1D motion along the reaction path and the orthogonal motions of the remaining (3N − 7) normal modes, is commonly referred to as reaction coordinates. It often corresponds to the dominant structural changes – in our case, the elongation and breaking of the O–H bond, while the remaining degrees of freedom correspond to vibrational modes orthogonal to this motion. The TS is located at s = 0, with negative and positive s values denoting the reactant and product sides, respectively.

At first, the 1D PES is constructed by computing the total energy V₀ at various points along the reaction path, while keeping the metal lattice rigid. Thereafter, the normal mode coordinates Q_k and corresponding frequencies ω_k(s) at those points are calculated via diagonalization of the force constant matrix. The eigenvectors, L_i,k(s), as obtained from the normal mode calculations, establish the connection between normal modes coordinates, Q_k and cartesian coordinates, x_i along reaction path coordinates s as:

where cartesian coordinates associated with the molecular configurations along the reaction coordinate are denoted by a_i(s).

The Hamiltonian corresponding to the reaction path coordinate can be written as:^88–90,93

wherewith the momenta conjugate to s and Q_k being represented as p_s and P_k, respectively. The other terms, b_ss and π_s are expressed as:

The B_k,j(s) term signifies the vibrationally non-adiabatic couplings with the following expression:

Note that the Hamiltonian is expanded to first order in b_ss and π_s, as higher-order terms are negligible. Energy exchange among the modes, k and j is governed by the Coriolis couplings terms, B_k,j through the π_s operator, while the interaction between a specific vibrational mode and the reaction coordinate (curvature coupling) is mediated by the b_ss operator via the B_k,9 terms. In the vibrationally adiabatic basis set, out of the nine normal modes for H₂O molecule, three of them correspond to the vibrational motion and the others to the rotational and translational motion. The “9th eigenvector” is particularly significant as it corresponds to the normalized gradient vector that characterizes movement along the reaction path. In other words, this eigenvector represents the mode with an imaginary frequency at the TS.

The total molecular wave function can be expanded in the close-coupled wave packet approach as the following:

where the adiabatic vibrational states of the molecule, Φ_v's are also the eigenfunctions of H_vib, with the corresponding eigenvalues being expressed by:

The index v in eqn (7) and (8) labels the vibrational state corresponding to the kth normal mode. The vibrationally adiabatic eigenfunctions, Φ_v's are product of 1D harmonic oscillator eigenfunctions having parametric dependence on s. By inserting the wave packets, Ψ(s;t) into the Hamiltonian given in eqn (2), the coupled equations-of-motion can be expressed in terms of χ_v(s;t) as:

The vibrationally adiabatic ground state (GS) is denoted by v = 0. In this present study, we consider only the single-quantum excitations together with the GS. The coupled wavepackets, in turn, propagate on the vibrationally adiabatic effective potential, V_eff, defined as:

These effective potentials simply denote the minimum energy paths (MEPs), where the zero point energy (ZPE) correction and the energy of any vibrational excitation along the reaction path have been properly incorporated. The operator F_vv' in eqn (9) consists of the b_ss and π_s operators which are responsible for the coupling of the vibrationally adiabatic states separated by single and multiple quanta. The strength of these couplings depend on the molecular momentum and kinetic energy. Consequently, at higher velocities and with larger coupling factors, curve crossings between vibrationally adiabatic potentials become more probable.

2.2 Propagation and analysis

The wavepacket is initialized in a specific vibrational state at the optimized geometry of the top-adsorbed H₂O molecule, where the vibrationally adiabatic PESs are almost constant with no coupling between them. A similar set-up with physisorbed H₂O had been used in previous theoretical studies^68,109 to motivate and reproduce quantum-resolved experiments. In our present calculations, the wavepacket is centered at a point in the reaction path with s ≈ −7.44 amu^1/2 Å in terms of mass weighted cartesian coordinates, a usual representation for the reaction coordinates, which translates to ≈−1.75 Å in the corresponding cartesian coordinates. The wavepacket is discretized on uniformly spaced grid having 512 grid points between −16.0 and 16.0 amu^1/2 Å. The effects of the kinetic energy and momentum operators on the wavepacket, χ_n are evaluated by the fast Fourier transform (FFT) technique.

We employ a cubic-form optical potential to absorb the wavepackets near the grid edges, preventing unphysical reflections from the boundaries, since the computational grid cannot extend to infinity in realistic quantum dynamical calculations. Such optical potentials, after being optimized, are applied at the last few grid points along the reaction path with |s| > 12.0 amu^1/2 Å and beyond the propagation time of ≈50 fs. The wavepacket is propagated with a time step of 0.01 fs until the total probability density on the grid falls below 10⁻⁶. In the product side, the non-adiabatic couplings become zero at s > 3.0 amu^1/2 Å and the adiabatic potentials are made constant beyond this point. Thereafter, the time-propagated wavepacket is Fourier transformed to extract energy-dependent amplitudes and consequently, energy resolved reaction probabilities, P^(k)₀(E_i,v) [= P^(k)_diss(E_i,v)] are computed from the reactive flux evaluated at large positive s values (5.5 and 7.0 amu^1/2 Å). Here, E_i denotes the set of incident energies included within the wavepacket, whereas v (= 0, 1 in present calculations) indicates the quantum number associated with the kth normal mode. The detailed mathematical frameworks underlying the evaluation of reaction probabilities are outlined in Section S1 of the SI. The dynamical parameters used in our calculations are tabulated in Table S1 of the SI.

3 Results & discussion

3.1 Normal mode analysis

We begin by examining the changes in the frequencies of the normal modes along the reaction path (s), as shown in Fig. 1, for H₂O dissociation on Pt(111) surface, both with and without 2D confinements. For all four systems representing different confinement scenarios, we observe that when the H₂O molecule is near the initial state, as indicated by negative s values, three of the normal modes, viz., symmetric stretching (ν_symm), asymmetric stretching (ν_asym), and bending (ν_bend), exhibit significant magnitudes of the frequencies (≈0.2–0.5 eV). Other five (5) modes (ν₄ to ν₈) correspond to the translational and rotational modes parallel to the surface, which become hindered at the TS and its vicinity on interaction with the surface. The frequencies of these modes show little variation, all remaining below 0.1 eV, indicating insignificant contributions of these modes to the reaction barrier. In the present study, we are primarily focused on the vibrational mode selectivity under different confining covers, and thus, the effects of those low frequency rotational and translational modes are not considered. Among the three vibrational modes, the frequency for the ν_asym mode remains nearly unchanged near the TS, for all the systems studied – bare Pt(111), Gr/Pt, C₃N₄/Pt and BN/Pt. This suggests that vibrational excitation of the ν_asym mode will have little effect on the reaction probability under adiabatic conditions. As a consequence, the reaction probabilities are expected to remain close to their respective GS values. In contrast, the frequencies of the ν_symm and ν_bend vibrational modes decrease remarkably at the TS. This behavior indicates that the excitation of these two modes could markedly enhance the reaction probabilities. In particular, the pronounced softening of the ν_symm mode within the vibrationally adiabatic approximation suggests an even greater potential for reactivity enhancement through its excitation. Notably, the extent of frequency softening in both the ν_symm and ν_bend modes is comparable across the different confined systems. However, in the case of C₃N₄ confinement, the softening for these two modes appears to occur more steeply along the reaction path between the reactant and the TS. Overall, the observed softening of vibrational frequencies underscores the profound impact of mode selectivity even within confined environments, reinforcing its role in steering reaction dynamics and energy redistribution through vibrational excitation.

Fig. 1 Normal mode frequencies of H₂O as a function of reaction path during H₂O dissociation on bare Pt(111), Gr/Pt(111), C₃N₄/Pt(111) and BN/Pt(111) systems.

3.2 Minimum energy path (MEP)

The GS potentials along MEP under different confined environments with and without zero-point energy (ZPE) correction are shown in Fig. 2. The ZPE-uncorrected barrier heights along the reaction path (s) for bare Pt, Gr/Pt, C₃N₄/Pt and BN/Pt systems are 0.92 eV, 0.69 eV, 0.81 eV and 0.80 eV, respectively. The obtained barrier and reaction energy values for bare Pt(111) (0.92 and 0.75 eV, respectively) are comparable to previous theoretical results [Table S3], taking into account methodological differences such as the choice of functional, surface coverage, and dispersion correction. This validates the reliability of the computed energy profiles in this work. Upon the inclusion of vibrational ZPE corrections for the respective systems, the barrier heights are diminished further by ≈0.14 eV [Table S2]. It is clearly evident in Fig. 3 that 2D confinement has a major role in reducing the barrier height for H₂O dissociation on Pt(111) surface with the most notable reduction observed for Gr, followed by BN. This observation can be attributed to the extent of electronic contributions from the 2D overlayers, as quantified by the Bader charge donation of 0.26e, 0.12e and 0.07e from Gr, h-BN and g-C₃N₄, respectively. The accumulation of the charge density in the antibonding molecular orbital of the H₂O molecule facilitates its dissociation.⁵⁹ The width of the MEP also decreases under confinement by the 2D layers. This trend suggests that the GS reactivity is likely to be increased for the confined cases. On pristine Pt(111) surface, the adsorption energy of H₂O is considerably smaller than the energy required for its dissociation;⁵⁵ hence, water splitting is not expected in the absence of surface promoters such as preadsorbed O or OH species. The present study, therefore, focuses on how two-dimensional confinement modifies the intrinsic O–H activation barrier once water is stabilized at the interface. The optimized structures of the TSs and other states corresponding to different s values along the minimum energy path for the dissociation of H₂O on Pt(111), Gr/Pt(111), C₃N₄/Pt(111) and BN/Pt(111) are also shown in Fig. 2. The distance between H₂O and the confining cover is about 2.7 Å. The length of the dissociating O–H bond in the initial state remains largely unaffected by confinement under different covers fixed at ≈5 Å above the metal surface [Table S4 in the SI]. This height corresponds to the minimum energy distance in presence of a H₂O adsorbate, with the attractive van der Waals forces and the repulsive interactions balancing each other out.⁶⁵ The same cover height has been used in the present calculations. In Fig. 4, the ZPE-corrected GS potential and the potentials corresponding to one-quantum excitation of ν_symm, ν_asym and ν_bend vibrational modes, are displayed. We collectively refer to these as the effective potential, V_eff [see eqn (10)]. As expected from the softening of the vibrational mode frequencies near the TS, the barrier heights associated with one-quantum excitation of ν_symm mode exhibit the highest decrease followed by those associated with the ν_bend mode [Fig. 4]. On the other hand, the profiles corresponding to one-quantum excitation of ν_asym mode are largely similar to those of ZPE-corrected GS potential, owing to its negligible softening near the TS as shown in Fig. 1. The introduction of confinement leads to a reduction in the barrier heights for all the specific mode-excited cases, with the extent of reduction following the order: Gr/Pt > BN/Pt ≈ C₃N₄/Pt > bare Pt. The ν_symm mode is most sensitive to confinement while the ν_asym mode is the least affected one [see Table S2]. Although the potentials for C₃N₄/Pt are much narrower in width, its barrier heights are only slightly higher in magnitude than those for BN/Pt. Therefore, the reactivity in the case of the C₃N₄ confinement is expected to resemble that of BN confined case.

Fig. 2 Comparison of ZPE-corrected and ZPE-uncorrected ground state potential for bare Pt(111), Gr/Pt(111), C₃N₄/Pt(111) and BN/Pt.

Fig. 3 Barrier heights for ZPE-corrected ground state potential and the potentials corresponding to one-quantum excitation of symmetric, asymmetric and bending vibrational modes.

Fig. 4 ZPE-corrected ground state potential and the potentials corresponding to one-quantum excitation of symmetric, asymmetric and bending vibrational modes for bare Pt(111), Gr/Pt(111), C₃N₄/Pt(111) and BN/Pt systems.

3.3 Reaction probabilities: adiabatic calculation

After the time propagation of the wavepacket under the influence of the respective effective potentials, we obtain the reaction probabilities for the concerned dissociation process of H₂O. Since our primary focus is to examine the influence of 2D confinement on the dissociative chemisorption process, the wavepacket is initialized between the Pt(111) surface and the 2D overlayer, where the interaction between the incoming H₂O molecule and the Pt(111) surface is minimized and simultaneously, the 2D confining layer warrants to have significant effect in modulating the effective potential profiles, as evident in Fig. 2 and 3. Also, it is noticeable from Fig. 2 and 4 that the magnitudes of the effective potentials are negligibly small beyond the |s|-value of ∼5.0 amu^1/2 Å for bare Pt and ∼3.0 amu^1/2 Å for the confined systems. Accordingly, the wavepacket for the present calculations are constructed by keeping it initially centered at sufficiently large negative s-value (≈−7.44 amu^1/2 Å) where the interactions with the V_eff are minimal and naturally, the couplings between the adiabatic PESs for various normal modes within a particular system are essentially close to zero. Importantly, systematic scanning of the center of the wavepacket confirm that initializing at even greater distance from the surface does not alter either the magnitudes or the trends of the calculated reaction probabilities, which clearly indicates the convergence of such profiles [see Fig. S2 and S3 and the associated discussion in Section S2 of the SI]. Fig. 5 presents the static surface dissociation probabilities (P_diss) of H₂O on Pt(111) with various 2D confinements, computed under the adiabatic approximation for the GS and for vibrationally one-quantum excited ν_symm, ν_asym and ν_bend modes. The bare Pt(111) surface shows the lowest reactivity for the GS as well as for the vibrationally one-quantum excited normal modes. Interestingly, the reaction probability profiles exhibit a sudden drop once the E_i decreases from the ZPE-corrected effective potential barrier; for example, in the GS case on bare Pt, the probability declines abruptly below ≈0.8 eV. This occurs because at energies lower than this threshold, the E_i is no longer sufficient to surmount the V_eff, and the reaction can proceed only via quantum tunneling. For the confined cases, on the other hand, the reactivities increase quite substantially with respect to (w.r.t.) the bare Pt(111) case. Under confinement by the Gr layer, the reactivity threshold decreases by more than ≈0.2 eV for both the GS and the mode-excited states, compared to the bare Pt(111) case. Just as observed for bare Pt(111), a sharp fall in the reaction probability appears beneath E_i corresponding to the ZPE-corrected effective potentials associated with the GS and one-quantum excited normal modes. In the case of C₃N₄/Pt(111), the threshold E_i, as required to reach reaction probabilities in the order of ≈10⁻⁴, diminishes by around ≈0.1 eV from the corresponding threshold for bare Pt(111). Notably, C₃N₄ confinement yields slightly higher reactivity than BN confinement for the GS as well as the one-quantum excited ν_asym and ν_bend modes [see Fig. 5], despite both overlayers resulting in comparable barrier heights, as tabulated in the Table S2. The narrower width of the MEPs in case of C₃N₄/Pt system likely contributes to its higher reactivity compared to the BN confined case, by facilitating enhanced quantum mechanical tunneling. Moreover, among all the systems studied, C₃N₄ confinement exhibits the narrowest MEPs, leading to a more gradual decline in P_diss and less steep profiles. A comparison of the mode-specific reactivities relative to the GS reveals that the one-quantum excitation of the ν_symm mode produces the highest reactivity, for cases with and without confinement, as shown in Fig. S4. This is due to the lower barrier height arising from the softening of the ν_symm frequency near the TS [Fig. 1]. In contrast, one-quantum excited ν_asym mode exhibits reaction probabilities similar to the GS, as it does not undergo any appreciable softening near s = 0. Therefore, within the static surface vibrational adiabatic condition, the reactivity order of various modes for a particular system are: ν_symm ≫ ν_bend > ν_asym ≈ GS [Fig. S4], whereas in terms of reactivity across different confined systems for a given one-quantum excited normal mode or the GS, the sequence is: Gr/Pt(111) > C₃N₄/Pt(111) ≈ BN/Pt(111) > bare Pt(111) [Fig. 5]. When the reactivity profiles of the ν_asym mode for the H₂O dissociation on bare Pt(111) surface is compared with that of systems involving other bare surfaces, viz., Ni,⁹³ or Cu and Ni–Cu alloy surfaces,⁷ a similar trend is observed. Remarkably, the reactivity achieved by exciting the ν_symm mode on bare Pt(111) (threshold ≈ 0.5 eV) is already attained at the GS under Gr confinement, demonstrating that confinement can substitute for mode-specific vibrational activation. Similarly, the reactivity corresponding to the one-quantum excited ν_bend mode on bare Pt(111) is closely reproduced at the GS in the BN-confined system. These observations highlight that the probability of H₂O dissociation on Pt(111) can be efficiently modulated through selective 2D confinement, thereby obviating the need for vibrational excitation. Although direct experimental comparison is not yet feasible for the 2D-confined cases, these results offer a new perspective on the dynamics of heterogeneous catalytic processes under confinement.

Fig. 5 Comparison of reaction probabilities for H₂O dissociation on different confined Pt(111)-systems as a function of incident energy of H₂O in the ground state as well as one-quantum excited vibrational modes within the vibrational adiabatic limit.

The vibrational efficacies,^73,89η^(k), of the three normal modes of H₂O, calculated for reaction probabilities up to ≈0.25, are displayed in Fig. 6. These profiles reflect the influence of different 2D overlayers on the mode-specific reactivity of Pt(111) catalyzed H₂O dissociation. For different normal modes, the η^(k) term is calculated as:

Fig. 6 Effect of 2D confinement on the vibrational efficacies for H₂O dissociation on Pt(111) for one-quantum excited (a) asymmetric stretching, (b) symmetric stretching, and (c) bending modes, evaluated within the vibrational adiabatic limit.

The expression in the denominator indicates the energy needed to excite the GS to the vth state of a particular mode k. In the numerator, E^(k)_i(v,P^(k)₀) and E^(k)_i(0,P^(k)₀) represent the incident translational energies required to achieve a specific dissociation probability P^(k)₀ for the vibrational state v and the GS (v = 0), respectively, for a particular vibrational normal mode k.

For the bare Pt(111) system, supplying an energy of ≈0.44 eV to the GS of the H₂O molecule facilitates a transition to the first excited state of the ν_asym mode [Table S5]. As shown in Fig. 6(a), the corresponding η values exhibit a small peak at low P_diss, before leveling off. Despite this non-monotonic trend, η remains consistently low (on average ≈0.02), indicating that vibrational excitation of the ν_asym mode has only a minor effect on improving the reaction probabilities. On the other hand, excitation of the molecule by ≈0.43 eV enables a one-quantum transition in the ν_symm mode. This results in remarkable enhancement in η, which rises to 0.50–0.53 for bare Pt(111) [Fig. 6(b)]. Similarly, the excitation of the GS with ≈0.19 eV, as necessary to reach the v_bend = 1 state, results in η ranging from 0.46 to 0.48 [Fig. 6(c)]. The η-profile for the ν_bend mode resembles that of the ν_asym mode, but with appreciably higher values. These observations clearly indicate that, within the vibrational adiabatic limit, excitation of the ν_symm or ν_bend mode has greater contribution in enhancement of reactivity, unlike the ν_asym mode. It is also emphasized that the relative efficacy of normal modes can be qualitatively inferred from the extent of their softening and the subsequent reduction in the activation barrier. Interestingly, under confinement by Gr and C₃N₄ overlayers, vibrational excitation of the ν_asym mode yields negative η values, implying a suppression of reactivity relative to the vibrational GS. Across the confined systems studied, η for the ν_asym mode remains nearly invariant with reaction probability and does not exceed 0.01. In the case of the ν_symm mode, the behavior of η varies considerably depending on the confining layer. The C₃N₄/Pt system shows the lowest η values up to a reaction probability of ≈0.18, beyond which it surpasses those of bare Pt. Confinement by Gr produces reasonably high magnitudes of η in the range of ≈0.64–0.66, followed by slightly lower values (≈0.60–0.63) observed for the BN/Pt system. For the one-quantum excited ν_bend mode, BN confinement maintains η values comparable to those of bare Pt, which exhibits the overall highest efficacy. The Gr overlayer leads to distinctly reduced η values (≈0.36–0.38) in comparison to other systems, suggesting a confinement-induced attenuation of vibrational promotion. The η values for C₃N₄/Pt (≈0.41–0.42) lie between those of BN/Pt and Gr/Pt, indicating a moderate contribution of its ν_bend mode in promoting reactivity. Notably, η appears relatively insensitive to P_diss beyond the initial low-dissipation regime, for most cases. In the presence of confinement, only the ν_symm mode displays enhancement in its ability to drive the process forward.

An intriguing aspect of this reaction is the emergence of nonstatistical behavior w.r.t. the E_i of the H₂O molecule. The behavior would be purely statistical if η were observed to be ≈1.0 across all the modes, which is evidently not the case. The obtained η values are less than 1.0, suggesting that vibrational excitation contributes less effectively to improving reactivity than translational energy. Despite the relatively limited contribution of vibrational modes, the strategic integration of mode selective excitation and 2D confinement offers a promising strategy for achieving precise control over catalytic H₂O dissociation.

3.4 Reaction probabilities: full coupling

Under experimental conditions, vibrational modes may couple with one another, enabling the redistribution of excess vibrational energy into translational motion along the reaction path, thereby allowing dissociation of H₂O molecule even at low E_i. To account for this energy transfer mechanism, it is essential to go beyond the vibrational adiabatic approximation and address the non-adiabatic effects by incorporating the Coriolis and the curvature coupling factors, which are represented by the mode–mode coupling (B_k,j) and mode–path coupling (B_k,9) terms, respectively [see eqn (4)–(6)]. Fig. S5 presents the reaction probabilities, obtained by including nonadiabaticity, for H₂O dissociation on various confined Pt(111) systems. The plots compare the probabilities for the GS and one-quantum excited vibrational modes. A key difference from the vibrationally adiabatic results is the markedly enhanced reactivity at E_i well below the corresponding ZPE-corrected barrier heights given in Table S2. Unlike the sharp decline seen in Fig. 5, the profiles from full coupling calculations exhibit a more gradual decrease with decreasing E_i, along with persistent zig-zag structures across all vibrational states [Fig. S5]. This behavior can be interpreted in terms of enhanced tunneling probabilities. The curvature of the reaction path allows the nuclear motion to redistribute into individual degrees of freedom. Such delocalization reduces the effective mass associated with the collective motion of the atoms along the reaction coordinate, thereby increasing the tunneling probability, particularly when curvature coupling is included.¹¹⁰

For the bare Pt(111) system, as illustrated in Fig. S5, the reaction probabilities for the one-quantum excited ν_asym mode approach unity above an energy of ≈0.78 eV. This value corresponds to the ZPE-corrected activation energy, as shown in Fig. 4, above which the reaction becomes classically allowed. Below this energy, if the molecule remains in the same vibrationally excited state (v_asym = 1), the reaction probability is expected to decrease due to insufficient translational energy to overcome the barrier. However, it is worth noting that, within a measurable likelihood, a molecule in a particular higher-energy vibrational state can undergo transitions to lower-energy vibrationally excited states or to the GS [see Fig. S6(b)]. These transitions can enhance reaction probabilities by converting excess vibrational energy into translational motion along the effective potential associated with the lower-energy modes [see Section S3 of SI].⁸⁶ Depending upon the coupling strength, the first excited state of the ν_asym mode can de-excite either to the first excited ν_symm or ν_bend mode, or directly to the GS. The ν_symm mode can relax either to the ν_bend mode or to the GS, whereas the ν_bend mode can make transition solely to the GS. These vibrational mode-to-mode and mode-to-GS de-excitation offer alternative pathway(s) that can promote reactivity, even when direct activation over the primary barrier is energetically unfavorable. We henceforth refer to these alternative pathway(s) as over-the-barrier (OTB) processes. The inclusion of nonadiabaticity increases the reaction probabilities by few orders of magnitude, particularly at lower E_i. Under fully coupled conditions, when both non-adiabatic coupling terms are included, threshold energy reductions relative to adiabatic profiles follow the trend: Gr/Pt(111) > BN/Pt(111) > C₃N₄/Pt(111) > bare Pt(111), corroborating the role of Gr overlayer in enhancing surface reactivity. The interplay of multiple de-excitation processes, classically allowed pathways, OTB processes, and quantum mechanical tunneling phenomena, along with their subtle variations, leads to the appearance of intricate features in the reaction probability profiles for all the normal modes, the details of which are discussed in the subsequent sections.

As shown in Fig. S5, the ν_symm mode with one-quantum excitation produces the highest P_diss for all the cases with or without confinement, similar to the adiabatic results. However, the inclusion of non-adiabatic effects leads to nearly barrierless behavior, which appears more prominently in the Gr-confined system, where the P_diss remains saturated over almost the entire energy range. This underscores the significant enhancement in vibrational reactivity brought about by full coupling. Another interesting outcome of including mode–mode and mode–path couplings is the appearance of a crossover between the reaction probability profiles of the ν_bend and ν_asym stretching modes, for all systems except Gr/Pt. At low E_i, the ν_asym mode is more reactive, but beyond a system-dependent crossover energy (E_cross), the ν_bend mode becomes dominant. The highest E_cross is observed for bare Pt(111). Accordingly, the general reactivity trend of the vibrational modes can be summarized as: (a) below E_cross:- ν_symm > ν_asym > ν_bend > GS; (b) above E_cross:- ν_symm > ν_bend > ν_asym > GS. These trends, obtained with the inclusion of non-adiabatic couplings, differ from the earlier findings on H₂O dissociation over Ni(111) as well as Cu(111) surfaces,^7,93 that had reported comparable reactivity for the one-quantum excited ν_asym and ν_symm modes when non-adiabatic terms were considered, despite the lower dissociation probability for one-quantum excited ν_asym mode under purely adiabatic conditions. The reactivity trend observed in the present case of bare Pt(111) surface can be attributed to the stronger value of ν_symm–ν_bend coupling, unlike the previous cases, where the ν_asym–ν_symm coupling possessed higher magnitude. Notably, for Gr/Pt, the reactivity corresponding to the vibrationally excited ν_bend mode also approaches a near-barrierless character and closely matches the profile of the ν_symm mode across the entire energy range. This distinct behavior, observed in the full-coupling calculations, is in stark contrast to the adiabatic results and other confined systems, where the ν_bend mode remains markedly less reactive than the ν_symm.

To unravel the impact of confinement, we examine how the reactivity associated with individual vibrational modes varies across different systems. Fig. 7, an alternate representation of Fig. S5, shows that the reactivities of the GS and the one-quantum excited ν_bend mode follow the ascending trend: bare Pt(111) < C₃N₄/Pt(111) ≈ BN/Pt(111) < Gr/Pt(111) in the low to intermediate energy regime (up to 0.6–0.7 eV), similar to the adiabatic results, whereas in the higher energy regime, the order for bare Pt(111) and C₃N₄ confined case gets interchanged. The ν_symm and ν_asym modes display deviation from this reactivity pattern, particularly at lower E_i. For the ν_symm mode, the reactivity order is: bare Pt(111) ≈ C₃N₄/Pt(111) ≈ BN/Pt(111) < Gr/Pt(111), with C₃N₄/Pt becoming the least reactive at higher E_i (0.6–0.9 eV). For the ν_asym mode, the initial trend is more atypical, with the reactivity order: Gr/Pt(111) < bare Pt(111) < C₃N₄/Pt(111) < BN/Pt(111). Above E_i ≈ 0.5 eV, however, this order realigns with the adiabatic expectations as well as the sequence observed for the GS and ν_bend modes. It is clearly noticeable that for the vibrationally excited ν_asym mode, the P_diss profile for the C₃N₄ confined system differs in overall magnitude from that for the BN confined case, unlike the behavior observed for the other modes, where similar P_diss profiles are obtained for these two systems. Moreover, profound oscillations are found in the ν_asym and ν_bend mode profiles for the C₃N₄ confined scenario. These distinctions highlight variations in vibrational energy transfer under different 2D confined environments and prompt further investigation into the underlying coupling mechanisms.

Fig. 7 Comparison of reaction probabilities for H₂O dissociation on different confined Pt(111)-systems as a function of incident energy of H₂O with the molecule being in the ground state or with one-quantum excited vibrational modes along with the incorporation of non-adiabatic couplings.

3.5 Non-adiabatic coupling

To delve deeper into the variations in reaction probabilities driven by non-adiabatic couplings, it is imperative to revisit the nature of the associated coupling elements. The B_k,9 term, responsible for reaction path curvature effects, plays a critical role in controlling the intramolecular energy flow, particularly in the vicinity of the TS. This curvature coupling facilitates interactions between the GS and singly excited vibrational states, thereby enabling efficient energy transfer between these states as the system evolves along the reaction coordinate. In parallel, the Coriolis coupling term, B_k,j, couples the singly excited states of the same symmetry. Thus, non-adiabatic coupling can increase the reactivity of the specific vibrationally excited modes by promoting the OTB dynamics, as discussed in the previous section. Within the close-coupled wave packet framework, both B_k,9 and B_k,j, derived from normal mode eigenvectors, are fitted to analytic forms and subsequently spline-interpolated onto the discrete uniform grid used for wavefunction propagation [Fig. 8].

Fig. 8 Curvature coupling values along the reaction path for the asymmetric stretching, symmetric stretching, and bending vibrational modes in (a) bare Pt(111), (b) Gr/Pt(111), (c) C₃N₄/Pt(111) and (d) BN/Pt(111) systems.

The subtle changes in B_k,9 as well as B_k,j values along the reaction path close to the TS, particularly before the TS, i.e., in the reactant region, are expected to be the major cause for the corresponding changes in the reactivity profiles. In Fig. 8, we have displayed the extent of mode–path couplings (B_k,9) along the reaction path for the ν_asym, ν_symm and ν_bend modes across four differently confined systems, viz., bare Pt(111), Gr/Pt(111), C₃N₄/Pt(111) and BN/Pt(111). It is clearly evident that the ν_symm mode exhibits significantly high magnitude of B_k,9 around the TS region for all the systems, whereas the other two modes generally possess negligible to moderate values, except for the Gr confinement, where the ν_bend mode is more strongly coupled to reaction path than the ν_symm mode. The couplings between different pairs of normal modes (B_k,j) as a function of the reaction path for the above four confined systems are shown in Fig. 9. Here, the magnitude for the coupling between the ν_symm and ν_bend modes appears to be much larger than the ν_asym–ν_bend and ν_asym–ν_symm coupling values.

Fig. 9 Coriolis coupling values along the reaction path for the asymmetric stretching, symmetric stretching, and bending vibrational modes in (a) bare Pt(111), (b) Gr/Pt(111), (c) C₃N₄/Pt(111) and (d) BN/Pt(111) systems. (The magnitudes along the ordinate have been rescaled by multiplication with a factor of ≈10³ for better visibility.)

For bare Pt(111), the ν_symm mode shows a strong coupling with the GS just before the TS, while the ν_asym and ν_bend modes, having negligible B_k,9 magnitudes, are weakly coupled with the reaction path [Fig. 8(a)]. Therefore, in case of non-adiabatic calculations with full coupling [Fig. 7 and Fig. S5], the one-quantum excited ν_symm mode shows a pronounced enhancement in reaction probability compared to the adiabatic profile [Fig. 5 and Fig. S4], thereby approaching barrierless behavior. Furthermore, energy transfer from the ν_symm to the ν_asym mode, facilitated by their coupling, contributes to an increase in P_diss for the latter. In contrast, one-quantum excitation of the ν_bend mode yields only a modest increment, although the ν_symm–ν_bend coupling is even stronger, probably because the only lower energy state accessible to ν_bend is the GS. For the Gr/Pt(111) system, the B_k,9 term is reasonably high for the ν_bend mode, in addition to the ν_symm mode [Fig. 8(b)]. This leads to a substantial increase in probabilities for the vibrationally excited ν_bend mode such that the probability profile of the ν_bend mode appears quite closer to the ν_symm mode, both depicting nearly barrierless nature [Fig. S5]. In the presence of C₃N₄ overlayer, the non-adiabatic reaction probabilities for the ν_symm mode are comparable to those on bare Pt(111) [Fig. 7], contrary to the relative increase predicted by adiabatic calculations [Fig. 5]. Meanwhile, the ν_bend mode for C₃N₄/Pt appears to be more reactive than bare Pt(111), particularly at low E_i, despite exhibiting negligible B_k,9 values like in case of bare Pt(111) [Fig. 8(c)]. Moreover, the relative ratio among different mode–mode couplings remain largely similar to those observed for bare Pt(111) although the confinement imposed by the C₃N₄ overlayer leads to marginally smaller magnitudes of the mode–mode coupling terms [Fig. 9(c)]. However, the relative decrease in barrier width and height along the MEP under C₃N₄ confinement w.r.t. bare Pt(111), as shown in Fig. 4 and Table S2, respectively, is slightly greater for ν_bend mode than for the ν_symm mode. Accordingly, the enhanced reactivity of ν_bend, along with the suppression in the expected increase in reaction probability for the ν_symm mode under C₃N₄ confinement compared to the bare Pt, can be attributed to the combined effects of variance in energetics and non-adiabatic couplings terms. In case of confinement by BN cover, along with the ν_symm mode, the ν_asym mode also exhibits noticeable coupling with the reaction path [Fig. 8(d)]. Also, notable is the relatively higher ratio of ν_asym–ν_bend coupling to ν_symm–ν_bend coupling in comparison with other systems. Hence, unlike the other confined systems, the reaction probabilities for the ν_asym mode in BN/Pt show greater amplification w.r.t. the bare Pt(111) case [Fig. S5], particularly in the lower energy regime, where the effects of non-adiabatic couplings are expected to have more influence nonetheless.

Thus, the enhanced reaction probabilities in the full coupling case can be rationalized from the concerted variations in B_k,9 and B_k,j values. These findings highlight the critical role of both mode–path and mode–mode couplings in governing mode selectivity during H₂O dissociation. Notably, the incorporation of 2D confinements can affect the strengths of curvature and Coriolis couplings, further influencing reactivity.

3.6 Individual effects of curvature coupling & Coriolis coupling

In Fig. 10, we isolate the individual effects of the non-adiabatic coupling parameters by examining two limiting cases in which each coupling term is selectively switched off. It is noticeable that when only the curvature coupling (B_k,9) is retained, in most cases the resulting trends and the structural features of the reaction probability profiles look quite similar to the ones obtained with full coupling calculations. In case of a singly excited ν_symm mode for all four systems, deactivation of the B_k,j term produces negligible changes in reaction probabilities, consistent with the dominant curvature coupling associated with this mode [see Fig. 8]. However, for the cases with considerable mode–mode coupling, the reactivities obtained with only the B_k,9 term are relatively smaller than those obtained from full coupling treatment. This effect is particularly evident in the one-quantum excited ν_asym and ν_bend modes for bare Pt(111), where significant deviations from the fully coupled profiles are observed upon neglecting B_k,j terms. These differences arise because of relatively stronger couplings of the ν_asym and ν_bend modes with the ν_symm mode, specifically in the reactant region of the reaction path, as can be seen in Fig. 9(a). Similar, though less pronounced, changes are observed in the reactivity profiles for the ν_asym and ν_bend modes for C₃N₄/Pt(111) system, due to lower values of ν_asym–ν_symm and ν_bend–ν_symm coupling strengths relative to those for bare Pt. Conversely, the reactivity for the systems confined by Gr and BN overlayers exhibit minimal sensitivity to the deactivation of the B_k,j terms, reflecting comparatively weaker mode–mode couplings. Therefore, when only curvature coupling is considered, the relative deviation in reaction probabilities, compared to those with full coupling can be arranged as: bare Pt(111) ≫ C₃N₄/Pt(111) ≈ BN/Pt(111) ≈ Gr/Pt(111). For both bare Pt(111) and C₃N₄/Pt(111), the order of variance in reactivity with various modes is as follows: ν_asym > ν_bend > ν_symm. The GS reactivity remains largely unaffected, as the curvature coupling enables energy flow to the one-quantum excited normal modes.

Fig. 10 Effect of different non-adiabatic coupling terms on the dissociation probabilities in the ground state as well as one-quantum excited vibrational modes on different confined Pt(111)-systems. “Solid lines” represent cases with only curvature coupling, whereas “dashed lines” denote cases with only Coriolis coupling.

A contrasting scenario emerges if solely the Coriolis coupling term is switched on: the reactivity trends resemble those predicted by vibrationally adiabatic calculations in cases where mode–path coupling dominates, but deviate from adiabatic behavior when mode–mode coupling becomes significant. For example, in all the studied systems, the highest increase in reactivity w.r.t. the adiabatic case [see Fig. S4 of SI] occurs in the ν_symm mode, even to the extent of being nearly barrierless. This can be explained from the fact that although switching off the substantially high B_k,9 coupling terms for the ν_symm mode is expected to result in similar probabilities like adiabatic ones, but mode–mode coupling involving ν_symm mode are comparatively larger than ν_asym–ν_bend coupling. Such an interplay between two kind of non-adiabatic coupling terms leads to the alteration in reactivity of the ν_symm mode to a lesser extent compared to the full coupling case. It is important to emphasize that both curvature and Coriolis couplings contribute to the emergence of additional classically allowed pathways, depending on their relative strengths. More specifically, mode–path couplings boost reactivity by making a transition to GS and thereafter, enabling the OTB processes and, also, initiating quantum mechanical tunneling processes. Introduction of mode–mode coupling, on the other hand, boosts reactivity to a lesser extent, whenever there is possibility of OTB processes via fetching a transition to energetically lower vibrationally excited states. Overall, the inclusion of only mode–path couplings result in much higher probabilities than inclusion of only mode–mode coupling, except in the case of ν_asym mode in bare Pt.

A careful examination of Fig. 10 reveals the complex characteristics of these extra classically allowed and quantum tunneling pathways in the bare Pt(111) system, arising on inclusion of only the mode–mode coupling terms. For visual clarity, the P_diss on the bare Pt(111) surface for various vibrationally excited modes, accounting solely for the Coriolis couplings, is presented in Fig. 11(a). The detailed structural features observed in the probability profile for the singly excited ν_asym mode, including plateaus and sharp declines, are interpreted in terms of OTB and tunneling processes, as annotated in Fig. 11(b). It is observed that in the case of singly excited ν_asym mode, the reaction probabilities reach unity for E_i above ≈0.80 eV, indicating efficient dissociation once the system can surmount the classical barrier height corresponding to v_asym = 1 [Table S2]. Below this energy, the reaction proceeds only via tunneling and there is a sudden drop in reactivity. Thereafter, a plateau region is observed for E_i ≈ 0.70–0.55 eV, indicating the onset of an OTB process involving v_asym = 1 → v_symm = 1 de-excitation (refer to Table S6). Below this energy, there is again a sharp decline in reactivity indicating the closure of the OTB pathway to v_symm = 1, with only viable mechanism being the tunneling through the barrier corresponding to this one-quantum excited ν_symm mode. In the energy range of 0.50–0.45 eV, a plateau region reappears, which can be identified as the OTB process via transition to the v_bend = 1 mode. For E_i in the range of 0.45–0.35 eV, again a sudden drop in reactivity is observed, which indicates that the OTB process through the ν_bend mode is closed. Next, for E_i < 0.35 eV, transition mainly to the GS occurs. For energy further below, reactivity is minimal due to closure of all classically allowed pathways with only feasibility being the quantum mechanical tunneling through the GS potential barrier. Likewise, for the one-quantum excited ν_symm mode, the first significant downturn in probability is observed for E_i < 0.55 eV, which corresponds to the incompetence to surpass the ZPE-corrected barrier height accompanying the one-quantum excited ν_symm mode [see Fig. 11(a)]. Subsequently at E_i ≈ 0.50–0.45, a plateau region emerges, reflecting the occurrence of OTB process by transition to the energetically lower one-quantum excited ν_bend mode. Thereafter, for E_i < 0.45 eV, again steep decline in reactivity is found signifying that the OTB process is barred to the one-quantum excited ν_bend mode, whereas for the energy range of ≈0.45–0.34 eV, the probabilities denote the occurrence of reaction via transition to the GS. The emergence of OTB processes via transitions from both v_asym = 1 and v_symm = 1 states to the GS occurs around similar E_i corroborating with the activation energy values in Table S6. Similarly, for the singly excited ν_bend mode, the probability profile approaches unity when E_i > 0.7 eV, which is the requisite energy to surmount the ZPE-corrected barrier height associated with this mode. Below this energy, there is sharp drop in reactivity as no other intermediate state is available for de-excitation except the GS [Fig. S6(b)], via which OTB process can occur for energies above 0.59 eV. Analogous physical insights emerge from the P_diss profiles of the confined systems, as depicted in Fig. 10. For instance, when the molecule resides in the first excited state of the ν_bend mode in case of Gr/Pt(111) system, the reaction probabilities tend toward unity beyond ≈0.56 eV as this energy surpasses the ZPE-corrected barrier. However, below this energy, the higher reactivity compared to the adiabatic results arises due to the lower activation energy (≈0.37 eV) for the OTB pathway, which becomes accessible following a transition to the GS. The profiles for the C₃N₄/Pt and BN/Pt systems also bear resemblance to the observations described above, with the OTB process initiating at a lower energy for BN/Pt than in C₃N₄/Pt [Table S6]. Hence, 2D confinement not only lowers the activation energy required for the classically permitted pathway on the vibrationally excited state, but also reduces the effective activation energy for OTB process following de-excitation to the energetically lower excited states or the GS.

Fig. 11 (a) Effect of inclusion of only mode–mode couplings on the reaction probabilities for H₂O dissociation on bare Pt(111) surface in the ground state as well as one-quantum excited vibrational modes as a function of incident energy. (b) For the one-quantum excited asymmetric stretching mode, depiction of alternative dissociation pathways in presence of only mode–mode coupling.

An interesting finding during this study is the relative prominence of these OTB features in the case of bare Pt(111) as compared to other confined environments, which can be attributed majorly to the relatively stronger mode–mode couplings in bare Pt. For the confined systems, mode–mode couplings are generally weaker, whereas the mode–path couplings are relatively stronger in comparison to those for bare Pt(111) system. Therefore, for the confined systems, switching off the mode–path couplings exhibit profiles like adiabatic results [Fig. 10]. Nevertheless, the OTB features are still visible for the one-quantum excitation of the ν_symm mode in the confined systems, except for Gr/Pt, where the corresponding activation energy is reduced to such an extent that these intricate OTB signatures disappear entirely. For the v_asym = 1 and v_bend = 1 states under confinement, the OTB features are much less prominent, particularly above P_diss ∼ 10⁻³. This can be rationalized by the comparatively weaker ν_asym–ν_symm and ν_asym–ν_bend couplings in the confined systems compared to those in bare Pt. In fact, when both mode–path and mode–mode couplings are included in the dynamical calculations [see Fig. 7 and Fig. S5], the OTB features get even more complicated with enhanced quantum mechanical tunneling due to lowering of effective mass in presence of curvature coupling effects. Thus, the bottom-line is that, without confinement, we can observe the detailed features of the OTB transitions to other energetically lower excited vibrational modes. However, for systems with confinements, owing to the weaker mode–mode couplings and strong mode–path couplings, transitions solely to the GS are observed below the energy corresponding to the ZPE-corrected barrier heights, a phenomenon that can be considered analogous to quenching to the GS.

4 Conclusion

In this study, we have thoroughly investigated the intricate dynamics of H₂O dissociation on Pt(111) surface, considering the influence of various two-dimensional (2D) confinements, by utilizing a comprehensive reaction path Hamiltonian approach. This work transcends from our previous theoretical investigation, where we had found that such 2D confinements using materials, viz., graphene, boron nitride, and graphitic carbon nitride induced different degrees of adsorption weakening and barrier reduction. Our present findings reveal the mechanistic subtleties behind vibrational mode-specific reactivity in both adiabatic and fully coupled (non-adiabatic) scenarios, especially under the influence of 2D confinements. The practical catalytic implications, though beyond the current scope, are likely to benefit from these mechanistic insights.

Through calculation of the reaction paths via construction of the minimum energy pathways, the incorporation of 2D overlayers appears to notably reshape the potential energy landscape by lowering the activation barriers, with Gr exhibiting the most pronounced effect. This observation is corroborated by the normal mode frequency profiles, where the symmetric stretching and bending modes exhibit substantial frequency softening as the system nears the transition state with the asymmetric mode remaining largely unaffected. Such frequency softening lowers the effective potential along the reaction path, highlighting the mode-specific reactivity trends by enabling easier energy redistribution and enhancing the likelihood of dissociation, particularly under quantum excitation. Furthermore, analysis of vibrational efficacies reveals that the symmetric stretching and bending modes have greater contribution towards the enhancement of reactivity relative to the asymmetric stretching mode.

Additionally, by incorporation of non-adiabatic effects, we gain a deeper understanding of the delicate interplay between curvature and Coriolis couplings, which govern the redistribution of vibrational energy into the reaction coordinate. These couplings not only amplify reactivity at low incident energies but also reveal non-statistical behaviors, as evident in the reactivity of the various vibrational modes. Through incorporation of 2D confinement, higher reaction probabilities are achieved, even to the extent where ground state reactivity is elevated to the levels comparable with those typically obtained for excited vibrational modes on bare surface. Among the 2D confined systems, Gr/Pt exhibits the most pronounced enhancement of probabilities, rendering the dissociation in case of singly excited symmetric stretching and bending modes nearly barrierless. This is attributed to strong mode–path and mode–mode couplings involving the symmetric and bending modes, with mode–path coupling having a greater impact on facilitating reactivity through extra classically-allowed pathways. Additional over-the-barrier processes arising from de-excitation to energetically lower excited states are also identified. The structural features originating from these processes are mostly prominent when mode–mode coupling dominates, as in the case of bare Pt, but are subdued for confined systems, where stronger curvature couplings overshadow the influence of Coriolis interactions.

To summarize, our findings highlight that 2D confinement allows for effective control over mode-specific dynamics of H₂O dissociation on Pt(111) surface. The insights gained here offer guidance for future experimental exploration and can be extended to other surface-mediated reactions, advancing our ability for rational design of catalytic systems where the tactical engineering of the overlayer architectures can be used to fine-tune vibrational energy flow and steer chemical reactivity with molecular-level precision. While the specific scenario of laser-prepared wave packets under confinement may be somewhat idealized, the qualitative mechanistic trends revealed by these calculations resonate with the previously observed effects of confinement.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available within the article and its supplementary information (SI). Supplementary information provides the details on numerical parameters for quantum dynamics, barrier heights for ZPE-corrected and -uncorrected potentials, evaluation of reaction probabilities, top and side view of the optimized initial states, and geometrical parameters of H₂O at different values of reaction coordinate. The convergence plots of the reaction probabilities under adiabatic condition are also included herein. The document also includes comparative plots of adiabatic and fully-coupled reaction probabilities across vibrational modes, estimation of accessible barrier heights from vibrational excitation, and tables summarizing the fundamental excitation energies for various vibrational modes of H₂O along with their corresponding over-the-barrier energy threshold. See DOI: https://doi.org/10.1039/d5cp03767k.

Acknowledgements

N. T. thanks Indian Institute of Science Education and Research (IISER) Kolkata for the Senior Research Fellowship. S. G. acknowledges ANRF-NPDF (File No.: PDF/2022/001824) as well as IISER Kolkata for financial grant. Both N. T. and S. G. are thankful to IISER Kolkata for the computational facility. A. K. T. acknowledges DST-ANRF New Delhi, India, for funding through project number CRG/2020/000040.

Notes and references

1. A. Vojvodic and J. K. Nørskov, Natl. Sci. Rev., 2015, 2, 140–143.
2. J.-H. Lin, P. Biswas, V. V. Guliants and S. Misture, Appl. Catal., A, 2010, 387, 87–94.
3. L.-Y. Gan, R.-Y. Tian, X.-B. Yang, H.-D. Lu and Y.-J. Zhao, J. Phys. Chem. C, 2011, 116, 745–752.
4. S. Ghosh, S. Hariharan and A. K. Tiwari, J. Phys. Chem. C, 2017, 121, 16351–16365.
5. M. García-Diéguez, I. Pieta, M. Herrera, M. Larrubia and L. Alemany, J. Catal., 2010, 270, 136–145.
6. S. Roy, S. Hariharan and A. K. Tiwari, J. Phys. Chem. C, 2018, 122, 10857–10870.
7. S. Ghosh, D. Ray and A. K. Tiwari, J. Chem. Phys., 2019, 150, 114702.
8. S. Roy, N. K. Jayadev, N. Tiwari and A. K. Tiwari, Int. Rev. Phys. Chem., 2020, 39, 267–318.
9. S. Roy and A. K. Tiwari, Phys. Chem. Chem. Phys., 2022, 24, 16596–16610.
10. X. Pan, Z. Fan, W. Chen, Y. Ding, H. Luo and X. Bao, Nat. Mater., 2007, 6, 507–511.
11. W. Chen, Z. Fan, X. Pan and X. Bao, J. Am. Chem. Soc., 2008, 130, 9414–9419.
12. M. Sanlés-Sobrido, M. Pérez-Lorenzo, B. Rodríguez-González, V. Salgueiriño and M. A. Correa-Duarte, Angew. Chem., Int. Ed., 2012, 51, 3877–3882.
13. J. Xiao, X. Pan, S. Guo, P. Ren and X. Bao, J. Am. Chem. Soc., 2014, 137, 477–482.
14. M. Zhao, K. Yuan, Y. Wang, G. Li, J. Guo, L. Gu, W. Hu, H. Zhao and Z. Tang, Nature, 2016, 539, 76–80.
15. G. Prieto, H. Tüysüz, N. Duyckaerts, J. Knossalla, G.-H. Wang and F. Schüth, Chem. Rev., 2016, 116, 14056–14119.
16. E. Fuentes-Fernandez, S. Jensen, K. Tan, S. Zuluaga, H. Wang, J. Li, T. Thonhauser and Y. Chabal, Appl. Sci., 2018, 8, 270.
17. X. Han, Q. Gao, Z. Yan, M. Ji, C. Long and H. Zhu, Nanoscale, 2021, 13, 1515–1528.
18. P. Sutter, J. T. Sadowski and E. A. Sutter, J. Am. Chem. Soc., 2010, 132, 8175–8179.
19. R. Mu, Q. Fu, L. Jin, L. Yu, G. Fang, D. Tan and X. Bao, Angew. Chem., Int. Ed., 2012, 51, 4856–4859.
20. L. Jin, Q. Fu, A. Dong, Y. Ning, Z. Wang, H. Bluhm and X. Bao, J. Phys. Chem. C, 2014, 118, 12391–12398.
21. M. Andersen, L. Hornekær and B. Hammer, Phys. Rev. B:Condens. Matter Mater. Phys., 2014, 90, 155428.
22. F. Yang, D. Deng, X. Pan, Q. Fu and X. Bao, Natl. Sci. Rev., 2015, 2, 183–201.
23. H. Li, J. Xiao, Q. Fu and X. Bao, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 5930–5934.
24. Q. Fu and X. Bao, Chem. Soc. Rev., 2017, 46, 1842–1874.
25. Q. Fu, W.-X. Li, Y. Yao, H. Liu, H.-Y. Su, D. Ma, X.-K. Gu, L. Chen, Z. Wang, H. Zhang, B. Wang and X. Bao, Science, 2010, 328, 1141–1144.
26. E. Grånäs, M. Andersen, M. A. Arman, T. Gerber, B. Hammer, J. Schnadt, J. N. Andersen, T. Michely and J. Knudsen, J. Phys. Chem. C, 2013, 117, 16438–16447.
27. Y. Yao, Q. Fu, Y. Y. Zhang, X. Weng, H. Li, M. Chen, L. Jin, A. Dong, R. Mu, P. Jiang, L. Liu, H. Bluhm, Z. Liu, S. B. Zhang and X. Bao, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 17023–17028.
28. M. Wei, Q. Fu, H. Wu, A. Dong and X. Bao, Top. Catal., 2016, 59, 543–549.
29. C. Zhao, J. Li, J. Dai, E. Voloshina, Y. Dedkov and Y. Cui, J. Phys. Chem. C, 2019, 123, 16137–16145.
30. I. Píš, E. Magnano, S. Nappini and F. Bondino, Chem. Sci., 2019, 10, 1857–1865.
31. S. Wang, Y. Feng, M. Yu, Q. Wan and S. Lin, ACS Appl. Mater. Interfaces, 2017, 9, 33267–33273.
32. L. Zhang, M. L. Ng and A. Vojvodic, J. Phys. Chem. Lett., 2020, 11, 9400–9407.
33. W. Zhao, B. Wang, M. Fan, L. Ling and R. Zhang, Mol. Catal., 2023, 548, 113421.
34. Y. Zhang, X. Weng, H. Li, H. Li, M. Wei, J. Xiao, Z. Liu, M. Chen, Q. Fu and X. Bao, Nano Lett., 2015, 15, 3616–3623.
35. L. Ferrighi, M. Datteo, G. Fazio and C. Di Valentin, J. Am. Chem. Soc., 2016, 138, 7365–7376.
36. Y. Zhou, W. Chen, P. Cui, J. Zeng, Z. Lin, E. Kaxiras and Z. Zhang, Nano Lett., 2016, 16, 6058–6063.
37. X. Chen and Z.-Z. Lin, J. Nanopart. Res., 2018, 20, 136.
38. X. Chen, Z.-Z. Lin, M. Ju and L.-X. Guo, Appl. Surf. Sci., 2019, 479, 685–692.
39. Z. Chen, C. Liu, X. Zhao, H. Yan, J. Li, P. Lyu, Y. Du, S. Xi, K. Chi, X. Chi, H. Xu, X. Li, W. Fu, K. Leng, S. J. Pennycook, S. Wang and K. P. Loh, Adv. Mater., 2019, 31, 1804763.
40. C. N. Eads, J. A. Boscoboinik, A. R. Head, A. Hunt, I. Waluyo, D. J. Stacchiola and S. A. Tenney, Angew. Chem., Int. Ed., 2021, 60, 10888–10894.
41. A. J. Shih, N. Arulmozhi and M. T. M. Koper, ACS Catal., 2021, 11, 10892–10901.
42. N. Tiwari, S. Hariharan and A. K. Tiwari, J. Chem. Phys., 2022, 157, 144701.
43. B. Wang, J. Diao, Q. Fu and Y. Ma, Appl. Surf. Sci., 2024, 656, 159715.
44. N. Tiwari, S. Ghosh and A. K. Tiwari, J. Phys. Chem. C, 2025, 129, 1988–1999.
45. M. Henderson, Surf. Sci. Rep., 2002, 46, 1–308.
46. C. Ratnasamy and J. P. Wagner, Catal. Rev., 2009, 51, 325–440.
47. C. Tsai, K. Lee, J. S. Yoo, X. Liu, H. Aljama, L. D. Chen, C. F. Dickens, T. S. Geisler, C. J. Guido, T. M. Joseph, C. S. Kirk, A. A. Latimer, B. Loong, R. J. McCarty, J. H. Montoya, L. Power, A. R. Singh, J. J. Willis, M. M. Winterkorn, M. Yuan, Z.-J. Zhao, J. Wilcox and J. K. Nørskov, Catal. Lett., 2016, 146, 718–724.
48. P. Ebrahimi, A. Kumar and M. Khraisheh, Emergent Mater., 2020, 3, 881–917.
49. W. T. Cahyanto, S. Zulaehah, W. Widanarto, F. Abdullatif, M. Effendi and H. Kasai, ACS Omega, 2021, 6, 10770–10775.
50. M. Mäkinen and K. Laasonen, Surf. Sci., 2023, 734, 122305.
51. L. Árnadóttir, E. M. Stuve and H. Jónsson, Surf. Sci., 2010, 604, 1978–1986.
52. M. J. Kolb, F. Calle-Vallejo, L. B. F. Juurlink and M. T. M. Koper, J. Chem. Phys., 2014, 140, 134708.
53. J. L. C. Fajín, M. N. D. S. Cordeiro and J. R. B. Gomes, J. Phys. Chem. A, 2014, 118, 5832–5840.
54. M. J. Ungerer, D. Santos-Carballal, A. Cadi-Essadek, C. G. C. E. van Sittert and N. H. de Leeuw, J. Phys. Chem. C, 2019, 123, 27465–27476.
55. G. B. Fisher and J. L. Gland, Surf. Sci., 1980, 94, 446–455.
56. X. Feng, S. Maier and M. Salmeron, J. Am. Chem. Soc., 2012, 134, 5662–5668.
57. Z. Li, S. Li, H.-Y. T. Chen, N. Gao, K. Schouteden, X. Qiang, J. Zhao, S. Brems, C. Huyghebaert and C. Van Haesendonck, J. Phys. Chem. Lett., 2019, 10, 3998–4002.
58. S. Ruiz-Barragan, D. Muñoz-Santiburcio and D. Marx, J. Phys. Chem. Lett., 2019, 10, 329–334.
59. F. Wei, Q. Wan, S. Lin and H. Guo, J. Phys. Chem. C, 2020, 124, 11564–11573.
60. M. J. Prieto, T. Mullan, M. Schlutow, D. M. Gottlob, L. C. Tănase, D. Menzel, J. Sauer, D. Usvyat, T. Schmidt and H.-J. Freund, J. Am. Chem. Soc., 2021, 143, 8780–8790.
61. M. Wang, C. Zhou, N. Akter, W. T. Tysoe, J. A. Boscoboinik and D. Lu, ACS Catal., 2020, 10, 6119–6128.
62. G. Xie, X. Liu, B. Guo, T. Tan and J. R. Gong, Adv. Mater., 2023, 36, 2211008.
63. E. Grånäs, U. A. Schröder, M. A. Arman, M. Andersen, T. Gerber, K. Schulte, J. N. Andersen, T. Michely, B. Hammer and J. Knudsen, J. Phys. Chem. C, 2022, 126, 4347–4354.
64. G. He, Q. Wang, H. K. Yu, D. Farías, Y. Liu and A. Politano, Nano Res., 2019, 12, 3101–3108.
65. N. Tiwari and A. K. Tiwari, ChemPhysChem, 2024, e202400586.
66. L. Juurlink, D. Killelea and A. Utz, Prog. Surf. Sci., 2009, 84, 69–134.
67. H. Chadwick, P. M. Hundt, M. E. van Reijzen, B. L. Yoder and R. D. Beck, J. Chem. Phys., 2014, 140, 034321.
68. P. M. Hundt, B. Jiang, M. E. van Reijzen, H. Guo and R. D. Beck, Science, 2014, 344, 504–507.
69. P. Kumar, H. Kuramochi, S. Takeuchi and T. Tahara, J. Phys. Chem. Lett., 2023, 14, 2845–2853.
70. B. Jiang, D. Xie and H. Guo, Chem. Sci., 2013, 4, 503–508.
71. B. Jiang and H. Guo, J. Chem. Phys., 2013, 139, 054112.
72. J. Li, B. Jiang and H. Guo, J. Chem. Phys., 2013, 139, 204103.
73. B. Jiang, M. Yang, D. Xie and H. Guo, Chem. Soc. Rev., 2016, 45, 3621–3640.
74. X. Zhou, F. Nattino, Y. Zhang, J. Chen, G.-J. Kroes, H. Guo and B. Jiang, Phys. Chem. Chem. Phys., 2017, 19, 30540–30550.
75. B. Jiang and H. Guo, J. Chem. Phys., 2014, 141, 034109.
76. B. J. Braams and J. M. Bowman, Int. Rev. Phys. Chem., 2009, 28, 577–606.
77. B. Jiang and H. Guo, J. Chem. Phys., 2015, 143, 164705.
78. T. Liu, Z. Zhang, B. Fu, X. Yang and D. H. Zhang, Phys. Chem. Chem. Phys., 2016, 18, 8537–8544.
79. T. Liu, Z. Zhang, B. Fu, X. Yang and D. H. Zhang, Chem. Sci., 2016, 7, 1840–1845.
80. Z. Zhang, T. Liu, B. Fu, X. Yang and D. H. Zhang, Nat. Commun., 2016, 7, 11953.
81. T. Liu, J. Chen, Z. Zhang, X. Shen, B. Fu and D. H. Zhang, J. Chem. Phys., 2018, 148, 144705.
82. T. Liu, B. Fu and D. H. Zhang, J. Chem. Phys., 2023, 158, 214305.
83. B. Jiang and H. Guo, J. Chem. Phys., 2013, 138, 234104.
84. B. Jiang, J. Li, D. Xie and H. Guo, J. Chem. Phys., 2013, 138, 044704.
85. B. Jiang and H. Guo, J. Phys. Chem. C, 2014, 118, 26851–26858.
86. B. Jackson and S. Nave, J. Chem. Phys., 2011, 135, 114701.
87. B. Jackson and S. Nave, J. Chem. Phys., 2013, 138, 174705.
88. M. Mastromatteo and B. Jackson, J. Chem. Phys., 2013, 139, 194701.
89. S. Nave, A. K. Tiwari and B. Jackson, J. Phys. Chem. A, 2014, 118, 9615–9631.
90. B. Jackson, F. Nattino and G.-J. Kroes, J. Chem. Phys., 2014, 141, 054102.
91. A. Farjamnia and B. Jackson, J. Chem. Phys., 2017, 146, 074704.
92. W. H. Miller, N. C. Handy and J. E. Adams, J. Chem. Phys., 1980, 72, 99–112.
93. H. Seenivasan, B. Jackson and A. K. Tiwari, J. Chem. Phys., 2017, 146, 074705.
94. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558–561.
95. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
96. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
97. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
98. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
99. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
100. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
101. S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem., 2011, 32, 1456–1465.
102. J. Carrasco, J. Klimeš and A. Michaelides, J. Chem. Phys., 2013, 138, 024708.
103. K. Tonigold and A. Groß, J. Comput. Chem., 2012, 33, 695–701.
104. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
105. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
106. T. Liu, T. Peng, B. Fu and D. H. Zhang, J. Phys. Chem. Lett., 2023, 14, 9713–9719.
107. T. Ohwaki, K. Kamegai and K. Yamashita, Bull. Chem. Soc. Jpn., 2002, 74, 1021–1029.
108. G. Henkelman, B. P. Uberuaga and H. Jónsson, J. Chem. Phys., 2000, 113, 9901–9904.
109. B. Jiang, X. Ren, D. Xie and H. Guo, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 10224–10227.
110. E. Kraka, Wiley Interdiscip. Rev.:Comput. Mol. Sci., 2011, 1, 531–556.

---

Footnotes

† These authors contributed equally to this work and share the first authorship.

‡ Present address: Department of Chemical Engineering, IIT Kanpur, Kanpur 208016, India.

§ Present address: Department of Chemistry, RRS College Campus, Patliputra University, Patna 803302, India.

---