Resolving sequential electron–proton transfer kinetics for electrochemical CO₂ reduction at the Cu(100)/H₂O interface via a quantum-classical framework

Abstract

Electron-transfer (ET) and proton-transfer (PT) events occurring at electrochemical interfaces ultimately dictate the efficiency of electrocatalytic energy conversion and chemical synthesis. Despite their significance, a unified theoretical description of ET/PT dynamics has been challenging due to the entangled electronic, nuclear, and solvent degrees of freedom. Herein, we present a quantum-classical multiscale framework integrating constrained density functional theory (CDFT) with machine learning accelerated molecular dynamics (MLMD) to investigate ET and PT during electrochemical CO₂ reduction on Cu(100) in explicit water. Distinct ML potentials are trained for the adiabatic ground state and for two charge-localized diabatic states, enabling efficient configurational sampling while preserving quantum-mechanical fidelity in electronic energies and forces. The derived diabatic free-energy surfaces reveal that CO₂ first undergoes inner-sphere ET to yield chemisorbed *CO₂⁻, followed by PT to form *COOH. Solvent reorganization imposes kinetic constraints on the ET step yet counterintuitively stabilizes the *CO₂⁻ intermediate through ion–dipole interactions, modulating the vibronic couplings. For PT, solvent relaxation dynamically adjusts the equilibrium donor–acceptor distance, thereby augmenting the Franck–Condon overlap between reactant and product vibronic wavefunctions in excited proton vibrational states, which facilitates nonadiabatic transitions across diabatic surfaces. Rate constants extracted by combining diabatic vibronic PCET theory with generalized Langevin equation-derived Grote–Hynes theory show that the sequential ET–PT pathway outpaces concerted PCET by about 5 orders of magnitude at the potential of zero charge (PZC). This methodology establishes a robust paradigm for dissecting ET and PT kinetics at electrochemical interfaces, emphasizing the interplay of quantum nuclear effects, vibronic couplings, and solvent fluctuations.

---

1 Introduction

The electrochemical reduction of CO₂ (CO₂RR) represents a compelling strategy for sustainably producing fuels and commodity chemicals, because it couples renewable electricity with CO₂ valorization and thus mitigates anthropogenic emissions.^1–5 Maximizing catalytic performance hinges on unequivocally identifying the rate-determining elementary step, which consensus attributes to coupled electron-transfer (ET) and proton-transfer (PT) events at the electrode–electrolyte interface.^6–9 Whether these transfers proceed sequentially via an inner-sphere ET followed by PT, or occur in a single concerted proton-coupled electron transfer (PCET) step, however, remains unresolved.^10–12 Conventional periodic density functional theory (DFT) calculations often favor a concerted PCET pathway for aqueous CO₂RR on metal/H₂O interfaces, whereas recent electro-kinetic studies point to a rate-limiting one-electron step whose kinetics are strongly modulated by the interfacial solvent structure.^13,14 This apparent discrepancy underscores the need for developing a theoretical framework that treats electrons, protons, and the dynamic solvent on an equal footing.

Several methodological limitations hinder such a framework. First, standard Kohn–Sham DFT is intrinsically adiabatic and cannot describe the nonadiabatic ET and PT processes that typify electrochemical reactions.^15–18 Constrained DFT (CDFT) can generate the requisite diabatic potential energy surfaces (PESs), yet exhaustive sampling of condensed-phase interfaces with CDFT remains computationally prohibitive. Second, thermodynamic integration (TI) proceeds along an alchemical coupling parameter rather than explicit collective variables, and thus ignores the stochastic coupling between the adsorbates, the polarizable metal surface, and the fluctuating hydrogen-bond network.^19–21 Third, kinetic models based on the Newns–Anderson Hamiltonian and Fermi's golden rule assume weak electronic coupling, an assumption that breaks down for chemisorbed intermediates that hybridize strongly with the Cu d-band.²² Finally, nuclear quantum effects (NQEs) – most notably proton tunnelling and zero-point motion can decisively influence kinetics but are routinely neglected in continuum or classical-MD solvent models.^23,24 Rectifying these deficiencies demands: (i) accurate diabatic landscapes that explicitly localize electronic charge; (ii) explicit, dynamical solvent models that capture reorganization energetics; and (iii) a quantum-mechanical treatment of proton motion to include tunnelling and kinetic isotope effects.

In this contribution, we deploy an integrated framework fusing CDFT with machine learning potentials (MLPs) to expedite molecular dynamics at the Cu(100)/H₂O interface. Distinct MLPs, trained for the adiabatic ground state and two charge-localized diabatic PESs, facilitate comprehensive configurational sampling while upholding quantum accuracy. Simulations reveal a prevailing sequential ET–PT mechanism: inner-sphere ET yields adsorbed *CO₂⁻, succeeded by PT to *COOH. A direct PCET pathway from physisorbed CO₂ is kinetically disfavored by a substantially higher barrier. Moreover, frequency-resolved analysis of the solvent friction kernel at the diabatic crossing uncovers specific interfacial vibrational modes that couple strongly to the reaction coordinate and dissipate energy. These findings demonstrate the necessity of multiscale, quantum-classical simulations that integrate electronic, nuclear, and solvent degrees of freedom to achieve a mechanistic understanding of electrochemical PCET processes.

2 Computational methods

2.1 Electronic structure

CDFT emerges as a robust methodology for extracting electron transfer (ET) parameters in condensed-phase systems, surpassing perturbative approaches by ensuring superior asymptotic behavior, global smoothness of diabatic potentials, and computational tractability.²⁵ In this work, CDFT-based molecular dynamics simulations of CO₂ at the Cu(100)/H₂O interface were executed within the CP2K software package.²⁶ The interface model comprised a 5 × 5 supercell with five Cu layers, described using the PBE-D3(BJ) functional.^27,28 The slabs were segregated by approximately 32 Å of explicit H₂O containing 170 molecules, maintaining an interfacial water density of about 1 ± 0.05 g cm⁻³ far from the surfaces.

2.2 Molecular dynamics

To efficiently explore configurational space and accurately capture interfacial energetics, deep neural networks (DNNs) were trained using energies and forces derived from CDFT calculations,^29,30 facilitated by a concurrent learning scheme within the ai2-kit toolkit.^31,32 Three separate MLPs were generated: one representing the adiabatic ground state and two corresponding to diabatic states localized with net charges of 0|e| and −1|e| on the CO₂, respectively. These charge-localized states were rigorously defined by optimizing constraint potentials using Lagrange multipliers within modified Kohn–Sham equations. The resulting MLPs effectively reproduced both CDFT and conventional DFT energies, ensuring accurate construction of diabatic and adiabatic free-energy surfaces (FESs). The reliability of these MLPs is demonstrated by their excellent performance on both training and test datasets, with detailed accuracy metrics provided in Fig. S2–S7. MD simulations employing these MLPs were carried out using the LAMMPS and PLUMED packages.^33,34 Simulations were conducted within the canonical ensemble (NVT), maintained at 330 K using a Nosé–Hoover thermostat, and integrated with a 0.5 fs timestep. Helmholtz free energies were calculated through well-tempered metadynamics (WT-MetaD), with Gaussian bias and collective variable parameters thoroughly detailed in the SI.

2.3 Rate theory

Accurate determination of PCET rate constants hinges upon the magnitude of vibronic coupling connecting reactant (µ) and product (ν) electron–proton vibronic states. These couplings are quantified by Hamiltonian matrix elements:³⁵

V_aµkν = V_akS_µν (1)

here V_ak is the electron coupling between reactant electronic state |a〉 and delocalized wavefunctions |k〉 of the quasi-free electrons in the metal, while S_µν represents the vibrational overlap integral between reactant and product proton vibrational states, i.e. 〈χµ^a|χ_ν^k〉. In the vibronically adiabatic regime, the PCET rate constant is accurately approximated by transition state theory (TST) modified by Grote–Hynes corrections^36,37 yielding:here, ΔG^‡ represents the free-energy barrier, k_B is the Boltzmann constant, T is the temperature, h is Planck's constant, and κ_GH is the Grote–Hynes transmission coefficient. However, within the diabatic vibronic formalism, the PCET rate explicitly depends on the vibronic coupling. For ET processes, the nonadiabatic rate constant can be expressed as:³⁸

For fixed proton donor-accept distance R, the Landau–Zener nonadiabatic PCET rate constants are evaluated as:^39–41

where w_µ is the Boltzmann probabilities for the proton vibrational state µ; is the probability of sampling a crossing point at . is the Maxwell distribution of velocities Ẋ, ω_X is the harmonic frequency along the collective reaction coordinate X, θ is the Heaviside step function, and is the cumulative weighted density of states for the occupied levels of the electrode. f(ε) is the Fermi–Dirac distribution function for the electronic states in the metal, ρ_M is the density of states at the Fermi level, β′ is a parameter representing the exponential decay of the electronic coupling with the distance between the molecule and the electrode, λ is the reorganization energy, and η is the applied potential. Notably, weak electronic coupling with macroscopic electrodes justifies the electronically nonadiabatic approximation per vibronic pair, permitting a truncated Landau–Zener expansion for transitions. Yet, the quasi-continuum of product states may confer apparent adiabaticity to the ensemble rate. Proton donor–acceptor motion is incorporated via thermal averaging:⁴²

k_tot = ∫dR P(R)k_CEPT(R) (5)

where P(R) are probability distribution functions for k_CEPT(R).

3 Results and discussions

3.1 Free-energy landscapes

The adiabatic and diabatic state free-energy landscapes of inner-sphere ET and CEPT are extracted from WT-MetaD simulations based on the MLPs. The primary reaction coordinate for inner-sphere ET is defined as the centroid distance along the z-axis between CO₂ and the nearest Cu(100) surface atom. For the subsequent protonation of *CO₂ to *COOH, two collective variables, R₁ and R₂ (illustrated in Fig. S1), are utilized to capture the proton transfer dynamics.

As shown in Fig. 1a, the FESs of adiabatic and diabatic state 2 (anionic CO₂⁻) exhibit nearly coincident minima in the product region. In contrast, in the reactant regime beyond a reaction coordinate of about 5 Å, the adiabatic surface aligns with diabatic state 1 (neutral CO₂). The essential distinction arises in the intermediate region, where the diabatic states intersect to form a crossing seam that mediates nonadiabatic transitions between electronic states. It is important to emphasize, however, that adiabatic and diabatic states are defined on distinct theoretical energy scales; therefore, direct comparison of absolute energies between the two is not physically meaningful. Crucially, the diabatic FESs intersect at 2.42 Å with a barrier of 0.37 eV, defining a nonadiabatic seam critical to interfacial ET dynamics. The displacement between the diabatic crossing (2.42 Å) and the adiabatic barrier (2.75 Å) reflects the additional free-energy cost associated with solvent reorganization required to stabilize the *CO₂⁻ species after the initial electronic interaction becomes favorable. Consistently, the net charge of CO₂ (Fig. S8) along the reaction coordinate also transition near the adiabatic barrier maximum. Following chemisorption and concurrent ET completion, PT from H₂O yields *COOH. The adiabatic ET–PT exhibits a reaction free energy of 0.27 eV and barrier of 0.34 eV. On the neutral diabatic surface-permitting only concerted PCET, the FES rises monotonically from reactants to products, devoid of a well-defined transition state and with a reaction free energy surpassing 1 eV, signifying that *COOH formation demands prior electron injection into CO₂. Conversely, the post-ET anionic diabatic surface displays *CO₂⁻ conversion to *COOH with a moderate barrier of 0.19 eV and free energy of −0.44 eV. Together, these observations substantiate a sequential ET–PT mechanism, wherein inner-sphere ET initially localizes charge on CO₂, enabling a feasible barrier for PT to yield *COOH. This framework clarifies how solvent electrostriction impedes the initial ET rate, while simultaneously accelerating protonation through ion–dipole stabilization of the anionic intermediate.

Fig. 1 Diabatic free-energy surfaces for CO₂ reduction on Cu(100). (a) One-dimensional profile for the inner-sphere ET step, plotted along the CO₂-surface centroid distance. (b–d) Two-dimensional free-energy landscapes for PCET projected onto the collective variables R₁ and R₂ on (b) the adiabatic ground state, (c) the neutral diabatic surface, and (d) the anionic diabatic surfaces.

3.2 Comparison of ET–PT and PCET kinetics

To elucidate kinetic distinctions between sequential ET–PT and concerted PCET pathways, we calculated the electrochemical ET and concerted PCET rate constants. Inner-sphere ET entails electron injection from Cu d-band states into the CO₂ lowest unoccupied molecular orbitals, i.e. 2p*. Electronic coupling at the diabatic crossing point was explicitly evaluated using the Newns–Anderson Hamiltonian (Fig. 2a), with the resulting values exceeding 0.38 eV. This contrasts with the lower-bound estimate of 0.11 eV obtained from a simplified Lorentzian fit (Fig. S11). This substantial coupling is consistent with strong chemisorption interactions. This large coupling indicates that the reaction proceeds on a single, ground-state adiabatic potential energy surface. In this regime, the reactant and the electrode share the transferring electron, the adiabatic TST is the appropriate framework for calculating the ET rate constant. While prior studies invoked Fermi's golden rule for ET rates in vacuo, solvent reorganization energy (λ) emerges as paramount at electrochemical interfaces.²² Structurally, forming CO₂⁻ entails a linear-to-bent geometric reconfiguration, with charge states eliciting differential solvent polarization.

Fig. 2 Electronic coupling and proton potentials at the crossing seam. (a) Electronic coupling matrix element V_ak obtained from the Newns–Anderson Hamiltonian fit at the centroid distance of 2.42 Å along the z-axis between the CO₂ molecule and the surface Cu atoms. (b) Proton potentials and corresponding vibrational wavefunctions for the contributing vibrational state for the reactant and product state at the equilibrium donor–acceptor distance R of 2.70 Å. The energy scale on the left corresponds to the proton potential, while the scale on the right corresponds to the vibrational energy levels of the proton.

Solvent reorganization energy for inner-sphere ET was quantified via dual MLP-driven molecular dynamics trajectories on neutral and anionic diabatic surfaces. Physisorbed (Q_R) and chemisorbed (Q_P) geometries are chosen as the initial configurations on the neutral and anionic surfaces, equilibrated for 10 ns with snapshots every 1 ps (10 000 frames each). The probability distributions of the vertical energy gap for each diabatic state, as shown in Fig. S9 and S10. The computed λ of 1.46 eV reflects substantial energetic costs from interfacial water reorganization, driven by ion-dipole interactions-key to charge stabilization. Marcus theory then estimates the crossing barrier as 0.47 eV from λ and reaction free energy. However, for the subsequent calculation of the ET rate constant, the free-energy barrier of 0.37 eV obtained directly from the WT-MetaD simulation was employed, as it provides a more complete and dynamical description of the reaction pathway. Extending this to PCET,⁴³ reactant geometries (proton on H₂O) and product geometries (proton on *CO₂⁻) were optimized with fixed other nuclei to yield λ_R = 1.76 eV and λ_P = 1.28 eV, respectively; the averaged λ ≈ 1.52 eV approximates the WT-MetaD reaction free energy of 0.27 eV. At the potential of zero charge (PZC, −0.5 V vs. SHE), k_ET ≈ 10⁶ s⁻¹-markedly subdued relative to vacuum estimates of 10¹⁵ s⁻¹.²² This disparity underscores the dominance of solvent reorganization, demonstrating that inner-sphere ET at metal/water interfaces cannot be captured by purely electronic descriptions but instead requires explicit treatment of bath-induced fluctuations and dissipative dynamics.

Within the quantum-mechanical framework of PCET,^18,44 the electron and proton are treated quantum mechanically, whereas the remaining nuclear coordinates are treated classically. The pivotal quantity dictating nonadiabatic PCET is the vibronic coupling V_aµkν, which, along with the proton vibrational wavefunctions, exhibits pronounced dependence on the proton donor–acceptor distance R. For each fixed R, one-dimensional diabatic proton potentials and associated vibrational wavefunctions are computed along the trajectory linking equilibrium proton positions in reactant and product states. The equilibrium R was ascertained via the potential of mean force (PMF) along this coordinate using umbrella sampling, yielding the normalized Boltzmann probability distribution (Fig. 3a) with an equilibrium value of 2.70 Å. It should be noted that the PMF increases along the distance R, but does not grow unbounded, thus, the Boltzmann factor continues to contribute a non-negligible probability across the sampled coordinate range. Subsequently, one-dimensional diabatic proton potentials were constructed at discrete values of R ranging from 2.50 to 3.10 Å in 0.1 Å increments. For each fixed R, molecular dynamics trajectories of 1 ns duration were propagated on the adiabatic ML-PES, employing the SHAKE algorithm to constrain R. Averaged structures obtained from these trajectories served as diabatic crossing-point geometries for subsequent CDFT calculations, providing a set of one-dimensional proton potentials connecting reactant and product equilibrium proton positions.⁴⁵

Fig. 3 Thermodynamic and kinetic for PCET. (a) Potential of mean force (PMF) along the proton donor–acceptor distance R (blue line) and the associated normalized Boltzmann probability distribution (red line). (b) Potential-dependent rate constants for inner-sphere ET (blue) and concerted PCET (red) obtained from Landau–Zener rate theory.

Fig. 2b illustrates diabatic proton potentials for both reduced and oxidized *CO₂ states at R = 2.70 Å, additional potentials and wavefunctions are provided in the Fig. S12. The proton potentials are notably asymmetric; the reduced state potential minimum is located near the oxygen atom of CO₂, whereas the oxidized state potential minimum is closer to the oxygen atom of H₂O. Analysis indicates that decreasing R distances enhance proton vibrational wavefunction delocalization, substantially increasing overlaps between reactant and product ground-state wavefunctions. Conversely, at larger R distances, excited vibrational states progressively dominate rate constant contributions due to increased overlaps with ground product states.⁴⁶

To accurately quantify the PCET rate constant, the individual R-dependent rate constants are integrated across the sampled range, each weighted by the Boltzmann probability distribution (Fig. 3a). At the equilibrium distance of 2.70 Å, significant contributions arise predominantly from transitions between the third excited reactant state and ground product state and from the first excited reactant and second excited product states. These excited-state transitions have substantial proton delocalization facilitating efficient vibronic coupling. Conversely, at a shorter distance of 2.50 Å, the primary contribution arises from ground-state transitions. Thus, rates at each R are governed by vibrational overlaps, favoring compact R; however, overall rates are modulated by sampling probabilities, biasing toward equilibrium R.⁴⁷ This interplay underscores the multifaceted balance shaping PCET kinetics.

Potential-dependent rate constants were evaluated using adiabatic TST for ET and the Landau–Zener rate expression for concerted PCET. The electrostatic potential at the reaction plane, relative to PZC was modeled via a reported Gouy–Chapman–Stern (GCS) approach (Fig. S13), incorporating distance-dependent dielectric permittivity in the outer Helmholtz layer to accurately represent interfacial dielectric saturation effects.⁴⁰ As illustrated in Fig. 3b, ET rate constants consistently exceed those of PCET across the examined potential range. This difference arises primarily from the inherently higher reorganization energy and diminished proton vibrational wavefunction overlaps intrinsic to PCET mechanisms. Kinetics are predominantly governed by ground-state vibrational transitions, capitalizing on substantial overlaps to enable efficient proton tunneling. Beyond the equilibrium R, excited vibronic states contribute approximately 25% to the PCET rate at the PZC when integrated over R distributions (Fig. S14). This mechanistic transition reflects a subtle balance among vibrational state populations, vibronic overlaps, and thermodynamic considerations. At shorter distances, excited states remain energetically prohibitive, leaving ground-state transitions driven by strong vibrational overlaps as the dominant contribution. At longer distances, despite the decreasing thermal accessibility and reduced conformational probability, excited vibronic states become increasingly favorable due to smaller free-energy gaps with the product state and substantial overlaps resulting from their vibrational delocalization.

Notably, the vibronic coupling strength initially increases with distance from 2.50 Å, driven by reduced vibrational energy level spacing and significant ground-state overlap integrals S₀₀. Beyond the equilibrium distance of 2.70 Å, the coupling strength diminishes, as contributions are increasingly dominated by delocalized excited vibronic states whose Boltzmann populations w_µ decline with increasing vibrational energy. Nevertheless, these excited-state transitions remain essential, contributing significantly (∼25%) to the overall PCET rate due to their enhanced thermodynamic favorability at larger distances. Thus, the fractional contribution of each vibrational state to the PCET rate emerges from a complex interplay between vibrational overlaps S_µν, the Boltzmann populations w_µ, and the state-specific free energy gaps ΔG_µ0.⁴⁵

3.3 Solvent dynamics

The utilization of equilibrium distribution functions W_µν(R) and M(Ẋ) in PCET theory for computing crossing-point sampling probabilities presupposes rapid solvent dynamics, wherein the solvent attains thermal equilibrium on a timescale much shorter than the reaction. This limit corresponds to the spatial diffusion regime of Kramers theory with Markovian friction. Nevertheless, polar solvent dynamics exert profound influence on reaction kinetics, particularly in processes entailing substantial charge redistribution at barrier crossings. Although curve-crossing formulas yield cusp-shaped barriers, solvent reorganization often emerges as the rate-determining factor. Furthermore, the Landau–Zener treatment presumes a constant velocity through the diabatic crossing; any velocity relaxation enhances adiabatic character and weakens nonadiabatic suppression. To mitigate solvent response delays stemming from inertial polarization, we adopt GH theory within the generalized Langevin equation (GLE) formalism, which embeds the reaction coordinate in a bath characterized by a frequency-dependent friction kernel. This framework is extended to PCET on diabatic FESs, with the friction kernel extracted from crossing point force autocorrelation functions (FACFs) computed from MLMD simulations, as depicted in Fig. 4. The friction kernel decays on a sub-100 fs timescale but exhibits pronounced oscillations, evidencing strong coupling between the reaction coordinate and interfacial solvent modes. The corresponding friction spectrum is dominated by a single band near 2500 cm⁻¹, attributable to the O–H stretching manifold of water. These ultrafast vibrations intermittently modulate the diabatic gap and shape the stochastic sampling of the crossing seam. The high-frequency intramolecular vibrations concurrently govern two critical aspects: they dictate the equilibrium sampling of the classical coordinate R by modulating the hydrogen-bond network, and they serve as the primary bath responsible for dissipating energy and imposing Markovian friction on the reaction dynamics. Incorporating the GH correction thus allows the approach to go beyond the instantaneous equilibrium (Markovian) assumption and reveals how inertial polarization dynamics dampen the overall charge transfer rate.

Fig. 4 Solvent dynamics at the PCET crossing seam. (a) Time-dependent friction kernel, governing the solvent reorganization along the reaction coordinate. (b) Corresponding friction spectrum, revealing dominant interfacial vibrational modes that couple to the reaction coordinate.

4 Conclusions

We have developed a comprehensive quantum-classical multiscale framework to elucidate ET and PT kinetics in electrochemical CO₂RR at the Cu(100)/H₂O interface, integrating CDFT, MLMD, and a diabatic vibronic rate theory augmented by GH corrections within the GLE formalism. The major conclusions are summarized as follows:

(1) Distinct neural network potentials for the adiabatic ground state and for two charge-localized diabatic states are well trained, allowing efficient configurational sampling while preserving quantum-mechanical fidelity for electronic energies and forces.

(2) Our calculations unambiguously discriminate between sequential ET–PT and concerted PCET pathways, showing that the sequential ET–PT pathway outpaces the concerted PCET pathway by approximately 5 orders of magnitude at PZC.

(3) Frequency-resolved friction analysis identifies interfacial O–H stretching modes as the principal bath coordinates that modulate electronic-nuclear coupling and the Landau–Zener transition probability.

(4) Vibrationally excited reactant states contribute 25% of the total PCET rate, highlighting the kinetic importance of excited-state vibronic channels in electrochemical environments.⁴⁸

(5) Acceleration of the rate-determining ET and PT events at the metal/aqueous interfaces can be achieved by stabilizing the *CO₂⁻ intermediate, by engineering hydrogen-bond networks that transiently shorten donor–acceptor separations, and by tailoring interfacial dielectric response to reduce solvent-induced friction.

Author contributions

Yun Yang: conceptualization, methodology, software, validation, formal analysis, investigation, data curation, writing – original draft, visualization. Gang Fu: conceptualization, resources, writing – review & editing, supervision, project administration, funding acquisition.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting the findings of this article are included in the supplementary information (SI). Supplementary Information: detailed information on the training errors of the MLPs, the CDFT setup, the electrical double layer model, the extrapolation of the vibronic coupling, friction kernel, and a comparison between the Newns–Anderson Hamiltonian and Lorentzian fitting for the electronic coupling matrix element, as well as additional computational details. See DOI: https://doi.org/10.1039/d5sc07385e.

Acknowledgements

We acknowledge support from National Key Research and Development Program of China (2022YFA1504601), the National Natural Science Foundation of China (22132004, U24A20489) and the Fundamental Research Funds for the Central Universities(20720250005).

Notes and references

1. M. Liu, Y. Pang, B. Zhang, P. D. Luna, O. Voznyy, J. Xu, X. Zheng, C. T. Dinh, F. Fan, C. Cao, F. P. G. de Arquer, T. S. Safaei, A. Mepham, A. Klinkova, E. Kumacheva, T. Filleter, D. Sinton, S. O. Kelley and E. H. Sargent, Nature, 2016, 537, 382–386.
2. Y. Wu, Z. Jiang, X. Lu, Y. Liang and H. Wang, Nature, 2019, 575, 639–642.
3. J. E. Huang, F. Li, A. Ozden, A. S. Rasouli, F. P. G. de Arquer, S. Liu, S. Zhang, M. Luo, X. Wang, Y. Lum, Y. Xu, K. Bertens, R. K. Miao, C.-T. Dinh, D. Sinton and E. H. Sargent, Science, 2021, 372, 1074–1078.
4. S. Ren, D. Jouli, D. Salvatore, K. Torbensen, M. Wang, M. Robert and C. P. Berlinguette, Science, 2019, 365, 367–369.
5. J. Gu, C. S. Hsu, L. Bai, H. M. Chen and X. Hu, Science, 2019, 364, 1091–1094.
6. I. V. Chernyshova, P. Somasundaran and S. Ponnurangam, Proc. Natl. Acad. Sci. U. S. A, 2018, 115, E9261–E9270.
7. A. Wuttig, M. Yaguchi, K. Motobayashi, M. Osawa and Y. Surendranath, Proc. Natl. Acad. Sci. U. S. A, 2016, 113, E4585–E4593.
8. W. Deng, P. Zhang, B. Seger and J. Gong, Nat. Commun., 2022, 13, 1–9.
9. J. A. Gauthier, M. Fields, M. Bajdich, L. D. Chen, R. B. Sandberg, K. Chan and J. K. Nørskov, J. Phys. Chem. C, 2019, 123, 29278–29283.
10. B. A. Zhang, C. Costentin and D. G. Nocera, J. Chem. Phys., 2020, 153, 094701–094707.
11. X. Qin, H. A. Hansen, K. Honkala and M. M. Melander, Nat. Commun., 2023, 14, 7607–7618.
12. M. Favaro, H. Xiao, T. Cheng, W. A. Goddard and E. J. Crumlin, Proc. Natl. Acad. Sci. U. S. A, 2017, 114, 6706–6711.
13. Y. Hori, in Electrochemical CO2 Reduction on Metal Electrodes, Springer, New York, 2012, vol. 488, pp. 89–189.
14. M. Gattrell, N. Gupta and A. Co, J. Electroanal. Chem., 2006, 594, 1–19.
15. A. A. Peterson, F. Abild-Pedersen, F. Studt, J. Rossmeisl and J. K. Nørskov, Energy Environ. Sci., 2010, 3, 1311–1315.
16. X. Nie, M. R. Esopi, M. J. Janik and A. Asthagiri, Angew. Chem., Int. Ed., 2013, 52, 2459–2462.
17. Y. C. Lam, A. V. Soudackov and S. Hammes-Schiffer, J. Phys. Chem. Lett., 2019, 10, 5312–5317.
18. C. Venkataraman, A. V. Soudackov and S. Hammes-Schiffer, J. Phys. Chem. C, 2008, 112, 12386–12397.
19. M. T. Koper, Chem. Sci., 2013, 4, 2710–2723.
20. A. J. Göttle and M. T. Koper, Chem. Sci., 2016, 8, 458–465.
21. X. H. Yang, A. Cuesta and J. Cheng, Phys. Chem. Chem. Phys., 2021, 23, 22035–22044.
22. S. Vijay, W. Ju, S. Brückner, S. C. Tsang, P. Strasser and K. Chan, Nat. Catal., 2021, 4, 1024–1031.
23. J. Waluk, J. Phys. Chem. Lett., 2024, 15, 598–607.
24. W. Zhu and Y. Zhao, J. Chem. Phys., 2007, 126, 184105–184113.
25. C. S. Ahart, K. M. Rosso and J. Blumberger, J. Chem. Theory Comput., 2022, 18, 4438–4446.
26. T. D. Kühne, M. Iannuzzi, M. D. Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Borštnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. Glöß, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack and J. Hutter, J. Chem. Phys., 2020, 152, 194103–194150.
27. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
28. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104–154123.
29. H. Wang, L. Zhang, J. Han and W. E, Comput. Phys. Commun., 2018, 228, 178–184.
30. J. Zeng, D. Zhang, D. Lu, P. Mo, Z. Li, Y. Chen, M. Rynik, L. Huang, Z. Li, S. Shi, Y. Wang, H. Ye, P. Tuo, J. Yang, Y. Ding, Y. Li, D. Tisi, Q. Zeng, H. Bao, Y. Xia, J. Huang, K. Muraoka, Y. Wang, J. Chang, F. Yuan, S. L. Bore, C. Cai, Y. Lin, B. Wang, J. Xu, J. X. Zhu, C. Luo, Y. Zhang, R. E. Goodall, W. Liang, A. K. Singh, S. Yao, J. Zhang, R. Wentzcovitch, J. Han, J. Liu, W. Jia, D. M. York, E. Weinan, R. Car, L. Zhang and H. Wang, J. Chem. Phys., 2023, 159, 054801–054825.
31. Y.-X. Guo, Y.-B. Zhuang, J. Shi and J. Cheng, J. Chem. Phys., 2023, 159, 094801–094813.
32. L. Zhang, H. Wang, R. Car and E. Weinan, Phys. Rev. Lett., 2021, 126, 236001–236008.
33. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171–108205.
34. G. A. Tribello, M. Bonomi, D. Branduardi, C. Camilloni and G. Bussi, Comput. Phys. Commun., 2014, 185, 604–613.
35. Y. Georgievskii and A. A. Stuchebrukhov, J. Chem. Phys., 2000, 113, 10438–10450.
36. R. F. Grote and J. T. Hynes, J. Chem. Phys., 1980, 73, 2715–2732.
37. M. Roca, V. Moliner, I. Tuñón and J. T. Hynes, J. Am. Chem. Soc., 2006, 128, 6186–6193.
38. E. Santos and W. Schmickler, Chem. Rev., 2022, 122, 10581–10598.
39. A. Soudackov and S. Hammes-Schiffer, J. Chem. Phys., 2000, 113, 2385–2396.
40. Y. C. Lam, A. V. Soudackov and S. Hammes-Schiffer, J. Phys. Chem. C, 2020, 124, 27309–27322.
41. R. E. Warburton, A. V. Soudackov and S. Hammes-Schiffer, Chem. Rev., 2022, 122, 10599–10650.
42. A. Soudackov, E. Hatcher and S. Hammes-Schiffer, J. Chem. Phys., 2005, 122, 014505–014518.
43. B. Auer, L. E. Fernandez and S. Hammes-Schiffer, J. Am. Chem. Soc., 2011, 133, 8282–8292.
44. P. Hutchison, K. Cui, J. Zhong and S. Hammes-Schiffer, J. Chem. Phys., 2025, 163, 091501–091516.
45. S. Horvath, L. E. Fernandez, A. V. Soudackov and S. Hammes-Schiffer, Proc. Natl. Acad. Sci. U. S. A, 2012, 109, 15663–15668.
46. S. Hammes-Schiffer, Energy Environ. Sci., 2012, 5, 7696–7703.
47. H. Ishikita, A. V. Soudackov and S. Hammes-Schiffer, J. Am. Chem. Soc., 2007, 129, 11146–11152.
48. B. J. Rousseau, X. Lu, A. V. Soudackov, J. Suntivich, H. Abruna and S. Hammes-Schiffer, J. Phys. Chem. C, 2024, 128, 12109–12120.

---