Computational advances in electrochemical CO₂ reduction: design of copper-based and atomically dispersed catalysts, interfacial dynamics, and machine learning

Conventional fixed-electrons DFT

Conventional fixed-electron DFT for eCO₂R is usually implemented through the CHE model^[16]. In this framework, one computes adsorption energies and intermediate stabilities on a slab with a fixed number of electrons, often in a simplified surface-vacuum setting, and then imposes the electrochemical driving force afterwards. In the original Nørskov formulation, the proton–electron pair is referenced to the standard hydrogen electrode (SHE) so that, at pH = 0 and 0 V vs. SHE, μ(H⁺ + e^-) is related to $$ \frac{1}{2} $$μ(H₂)^[16]. The effect of bias enters through a -eU_SHE term, with U_SHE the electrode potential referenced to the SHE, and pH is incorporated through the chemical potential of the proton. Relative to the pH = 0 reference, each proton involved in an elementary step contributes a correction of k_BTln10 × pH to the reaction free energy, corresponding to about 0.059 eV per pH unit at 298 K. As a result, proton-consuming steps become less favourable with increasing pH on the SHE scale, whereas this dependence is largely absorbed into the reference electrode when energies are expressed vs. the reversible hydrogen electrode (RHE), though this cancellation holds strictly only for concerted proton–electron transfer steps, not for chemical or electron-only steps whose energetics remain pH-independent on both scales. In this sense, CHE converts conventional surface DFT into potential- and pH-dependent reaction thermodynamics without explicitly simulating ions or a constant-potential electrode. CHE therefore makes it possible to build free-energy diagrams, identify the most uphill elementary step, estimate a theoretical limiting potential, U_L, defined as the least negative potential at which all electrochemical elementary steps become exergonic, U_L = -max_i(ΔG_i)/e, and map thermodynamic trends in adsorbate stability and product selectivity across catalyst families. These capabilities come with an important caveat: CHE cannot describe kinetic barriers, explicit solvent structure, or potential-dependent mechanism switching.

CHE is much less reliable when the question depends on explicit solvent structure, cation-specific effects, field-sensitive intermediates, decoupled sequential proton-transfer and electron-transfer steps, in which the two transfers do not occur concertedly, or quantitatively accurate activation barriers at a fixed electrode potential, or coverage-dependent lateral interactions between co-adsorbed intermediates whose competition determines surface populations under steady-state conditions; the implications of these limitations for eCO₂R mechanism studies are discussed further in Section “Beyond the CHE”^[17-19]. Therefore, in the context of eCO₂R, conventional fixed-electron DFT with CHE is best viewed as a first-pass thermodynamic framework for promising catalyst screening rather than a complete high-accuracy electrochemical theory to resolve detailed catalytic mechanisms at realistic electrochemical interfaces.

Grand canonical and constant potential DFT

While CHE-based fixed-electron DFT is useful for constructing approximate free-energy landscapes, its inability to maintain a specific electrode potential limits its accuracy for genuinely potential-dependent electrochemical interfaces. This motivates the use of grand canonical DFT (GC-DFT). This method treats the electrode as an open electronic system coupled to an external electron reservoir, representing the potentiostat in an electrochemical cell, such that the electron chemical potential μe = -eU_SHE is held fixed, where U_SHE is the applied electrode potential referenced to the SHE, and the number of electrons adjusts self-consistently to maintain the target electrode potential^[20]. In practice, GC-DFT fixes the electrode potential by holding the electron chemical potential μe = -eU constant while the number of electrons varies self-consistently; this constant-potential (grand canonical) treatment is the counterpart to a constant-charge (canonical) calculation, in which the electron number is fixed, and the potential is left to float. Realising a fixed potential requires the net charge of the cell to be neutralised and an absolute potential reference to be established, which is the role of the counter-charge model: the standard and more physical choice is an implicit Poisson–Boltzmann electrolyte, in which the compensating charge forms a diffuse ionic layer in a continuum electrolyte that also references the potential to the bulk solution, while a homogeneous (jellium) background charge is a cruder alternative that lacks such a reference. Implicit-electrolyte models are thus themselves a standard means of imposing constant potential, and the constant inner potential (CIP) formulation^[20], which constrains the electrostatic potential at the electrode rather than the Fermi level, is an equivalent constant-potential scheme rather than an alternative to the counter-charge representation. The methodologically important distinction is therefore between constant-charge and constant-potential treatments, not between implicit-electrolyte counter-charge and the constant-potential formulations. That makes it better aligned with operational conditions in electrochemistry than conventional fixed-electron DFT. Moreover, GC-DFT can directly determine how adsorbate geometries, surface charge, reaction kinetics, and energetics change with electrode potential rather than relying on a linear estimation as the CHE model does.

Many studies indicate that the non-linear response between electrode potential and intermediate geometries, binding strength, and charge distribution may be neglected by an excessively simplified CHE model and lead to misleading results^[21]. Non-linear dependence of intermediate geometry, binding strength, and surface charge on electrode potential has been documented in several eCO₂R-relevant systems^[22,23]. Simplifying this response to a linear -eU correction, as in CHE, can produce qualitatively incorrect predictions of the rate-determining step and preferred pathway. Furthermore, the CHE potential correction enters as a strictly linear -eU term applied uniformly to each electron-transfer step; it does not capture the non-linear dependence of intermediate geometry, binding strength, and surface charge on electrode potential that has been documented computationally for several eCO₂R-relevant systems. Fixed-electron DFT combined with transition-state search methods such as nudged elastic band (NEB) can provide first-order estimates of activation barriers at a given surface geometry, but the barrier is evaluated at a fixed electron count and the potential dependence is introduced only through a post hoc correction; however, this approach omits the non-linear response of the electric double layer and the effect of the potential of zero charge (PZC) on barrier heights, effects that are most self-consistently captured within a grand canonical framework. Additionally, the specific reaction mechanism may differ with the applied potential; for example, *CO dimerisation is more favourable at mild potentials while CO–CHO* coupling takes place at more reduced potentials for C₂ product formation^[24].

Grand canonical DFT therefore provides a more realistic description of electrochemical interfaces under applied bias than static fixed-electron approaches. Nevertheless, grand canonical DFT still treats the surface composition as fixed during each calculation; the fluctuating populations of co-adsorbed intermediates under steady-state electrochemical conditions are more naturally addressed by grand canonical Monte Carlo (GCMC) sampling, introduced in the following subsection.

GCMC

GCMC methods are an ensemble-sampling method in which the catalyst surface exchanges adsorbates with a reservoir through stochastic insertion and deletion moves, accepted or rejected according to a Metropolis criterion based on DFT-computed binding energies, so that the number and identity of adsorbed species fluctuate at fixed temperature and chemical potentials^[25]. The electrochemical driving force enters through the chemical potential of each electroactive adsorbate, which is coupled to the electrode potential U via μ_i(U) = μ_i⁰ - n_ieU_SHE, where U_SHE is the electrode potential referenced to the SHE, n_i is the number of electrons transferred in forming species i from solution-phase reactants; varying U therefore shifts the relative stability of all adsorbates simultaneously and yields a potential-dependent surface coverage directly comparable to operando spectroscopic measurements. In electrocatalysis, the primary value of GCMC is not the calculation of elementary reaction barriers, but the identification of realistic catalyst structures, surface coverages, and active-site populations under operating conditions^[25]. The dynamic competition between adsorbates can reveal active sites inaccessible to static DFT simulations^[26], explicitly capture potential-dependent co-adsorption and site-blocking^[27], and provide mechanistic insight into adsorbate-induced surface reconstruction^[28]. In eCO₂R studies, the GCMC-DFT integrated method has been used to investigate surface reconstruction under operational conditions and uncover hidden adsorption sites that can be missed^[27]. Additionally, large-scale GCMC simulations have been used to sample adsorbate-induced Cu surface roughening under eCO₂R conditions, linking high *H coverage to facet-dependent reconstruction and subsurface H incorporation^[29]. GCMC provides statistically sampled, condition-dependent structures and coverages that can feed into DFT energetics, micro-kinetics, and selectivity analysis^[30] [Figure 1]. The operando application of GCMC to Cu-based and single-atom catalyst (SAC) systems is discussed in detail in Section “Operando evolution and co-adsorption”. Nevertheless, GCMC samples configurational space stochastically and does not propagate nuclear motion in time; capturing the real-time dynamics of interfacial solvent, ion transport, and proton transfer requires MD-based approaches, introduced in the following subsection.

Figure 1. Coupled GCMC–EKMC workflow for simulating potential-dependent *CO/*H co-adsorption and Cu(100) surface reconstruction under eCO₂R conditions. Top left: Initial ordered Cu(100) surface (~35 nm lateral dimension)^[30]. Top right: The GCMC loop samples the occupancy of each surface site (empty, *CO-occupied, or *H-occupied) under fixed electrochemical parameters (U, pH, T, pCO) via Metropolis-criterion acceptance until convergence: a proposed move is accepted with probability min{1, exp(-ΔG_tot/k_BT)}, where ΔG_tot = G_totⁱ⁺¹ - G_totⁱ, so that moves lowering the total free energy are always accepted, and those raising it are accepted with Boltzmann probability^[30]. Bottom left: Final reconstructed Cu(100) surface (~35 nm), in which adsorbate-driven Cu atom migration produces a disordered morphology qualitatively distinct from the initial slab^[30]. Bottom right: The converged GCMC coverage is passed to an EKMC loop, which propagates the surface in time by distinguishing fast (CO and H diffusion) from slow (Cu migration) events^[30]. All panels reproduced with permission from Ref.^[30]. Copyright © 2026 American Chemical Society. GCMC–EKMC: Grand canonical Monte Carl and environmental kinetic Monte Carlo; eCO₂R: electrochemical CO₂ reduction.

Dynamic modelling

Even when realistic adsorbate coverages and surface compositions are established, for example through GCMC sampling, the electrochemical interface remains intrinsically dynamic: solvent molecules continuously reorient, hydrogen-bond networks fluctuate, and ions migrate on picosecond-to-nanosecond timescales that are inaccessible to static DFT or stochastic Monte Carlo methods. Capturing these effects requires MD simulations, in which nuclear positions are propagated in time according to forces derived either directly from DFT [ab initio molecular dynamics (AIMD)] or from a parameterised interatomic potential. In eCO₂R, AIMD is the most rigorous MD approach, providing explicit treatment of electronic structure, proton transfer, and bond breaking at the interface, but its high computational cost restricts accessible simulation times to tens of picoseconds and system sizes to a few hundred atoms^[17]. Key interfacial quantities such as water orientation distributions, hydrogen-bond lifetimes, ion pairing at the double layer, and proton transfer rate constants are obtained as time-averaged observables over these trajectories and cannot be inferred from single static configurations^[31]. Dynamic simulations therefore become essential when the chemical processes of interest depend on explicit solvent structure, ion-specific effects at the double layer, or kinetic bottlenecks involving proton or electron transfer steps that are sequential rather than concerted. However, the high cost of AIMD severely limits the length and timescales accessible, motivating the development of MLFFs that can reproduce AIMD-quality forces at a fraction of the computational cost, as described in the following subsection.

AIMD provides DFT-level accuracy by propagating nuclear motion at the solvated electrochemical interface, with trajectories typically spanning 10-100 ps over systems of a few hundred atoms^[32]. In eCO₂R, AIMD captures explicit water effects, proton shuttling through hydrogen-bond networks via the Grotthuss mechanism, dynamic adsorption and desorption of intermediates, short-range cation–intermediate interactions, and, when combined with enhanced sampling techniques such as metadynamics or umbrella sampling, free energy profiles and activation barriers along solvent-coupled reaction coordinates; all of these are inaccessible to static DFT calculations^[33]. AIMD therefore provides the most chemically realistic description of the eCO₂R interface currently accessible to first-principles simulation, at the cost of trajectories limited to tens of picoseconds and systems of a few hundred atoms^[34].

Given the computational cost of AIMD, trajectories are typically limited to 10-100 ps, often too short to directly observe key reaction steps such as intermediate protonation or C–C coupling, which occur on timescales inaccessible to unbiased simulation^[32]. To address this limitation, metadynamics (MetaD) provides an effective approach to accelerate rare events by reconstructing the free energy surface (FES) along selected collective variables through the accumulation of a history-dependent bias potential that progressively fills free energy minima and drives the system across otherwise rarely visited transition states^[35]. The key advantage of MetaD is that it transforms qualitative structural trajectories into quantitative free energy profiles and activation barriers for elementary steps too infrequent to be sampled within an unbiased trajectory. However, converged free energy profiles require long MetaD simulations that can themselves exceed the reach of AIMD, because convergence demands exhaustive sampling of all accessible free energy basins, a condition that becomes increasingly demanding as the complexity of the free energy landscape or the number of collective variables grows.

Classical force fields extend the spatial and temporal reach of atomistic simulation by replacing quantum mechanical force evaluation with analytical interatomic potential functions, eliminating the electronic structure calculation at each timestep and enabling trajectories of nanoseconds to microseconds over systems of millions of atoms^[36-38]. Standard non-reactive force fields, such as TIP4P for water or the embedded atom method (EAM) for metal surfaces, assume a fixed bonding topology and cannot describe bond breaking or formation^[39]. Reactive force fields such as ReaxFF overcome this restriction by using a bond-order formalism in which connectivity is determined dynamically from interatomic distances, allowing chemical reactions to be simulated at classical cost^[40-42]. In eCO₂R, ReaxFF has been applied to long-time equilibration of the electrode–electrolyte interface and to investigating how surface modifiers influence interfacial water organisation and hydride dynamics, problems that are computationally intractable with AIMD. However, the accuracy of ReaxFF for bond-breaking chemistry depends critically on the quality of the parametrisation, which is typically derived from gas-phase or vacuum quantum mechanical reference data and therefore does not capture the potential-dependent charge redistribution, surface polarisation, or solvation effects that govern reactivity at electrochemical interfaces^{[34,41,43,44]}. This transferability gap motivates the development of MLFFs, which retain the reactive flexibility and near-classical computational cost of empirical potentials while being trained directly on AIMD-quality data for the relevant interfacial environment, as described in the following subsection.

Machine learning forcefields

The spatial and temporal limitations of AIMD and enhanced-sampling approaches make it difficult to access the larger systems and longer trajectories needed to describe electrochemical interfaces under realistic operating conditions. MLFFs offer a practical route to extending such simulations by training interatomic potentials on DFT energies and forces, reproducing the potential energy surface with near-DFT accuracy at a fraction of the computational cost of direct AIMD^[45,46]. A key strategy for generating MLFFs is the on-the-fly active learning approach implemented in VASP by Jinnouchi et al., in which the force field is updated during the simulation whenever the Bayesian uncertainty of the predicted energy exceeds a threshold, triggering a new DFT single-point calculation and automatically expanding the training set to include the new configuration, making the approach well suited to reactive, fluctuating interfaces whose configurational space cannot be anticipated in advance^[47,48]. Beyond on-the-fly active learning, other widely used MLFF strategies include offline active-learning potentials^[49], deep-potential models^[50], equivariant graph neural network potentials such as MACE^[51], and in the electrochemical context, constant-potential MLFFs adapted to charged electrochemical interfaces^[52]. All of these retain near-DFT fidelity while extending accessible length and time scales far beyond conventional AIMD. In eCO₂R, MLFFs are particularly valuable when the catalyst–electrolyte interface is too dynamic for a static slab picture. Tian et al. applied an MLFF to the Cu–water interface under applied cathodic potential and found that potential-driven reorganisation of interfacial water and modification of the hydrogen-bond network measurably alter the barriers of key hydrogenation and ion-pair dissociation steps^[31], effects entirely inaccessible to static slab models and only marginally resolvable within the timescales of conventional AIMD. Compared with AIMD, MLFFs make it practical to simulate larger supercells, longer trajectories, explicit solvent environments, and more extensive rare-event sampling, capturing solvent reorganisation, double-layer fluctuations, potential-dependent interface dynamics, competitive C–C coupling pathways, and large adsorbate-configuration ensembles. A current limitation is that most MLFFs are trained at fixed electron count rather than fixed electrode potential, which restricts their fidelity for genuinely potential-dependent phenomena; the development of constant-potential MLFFs represents one of the most important near-term methodological priorities in the field^[52]. A further limitation shared by most current MLFFs is their use of a strictly local energy decomposition, in which each atomic contribution depends only on the chemical environment within a finite cutoff radius. This approximation neglects long-range electrostatic interactions from the charged double layer, restricts explicit representation of charge transfer between the metal surface and adsorbed intermediates, and cannot capture non-local through-metal coupling between co-adsorbed species, all of which contribute to adsorption energies and lateral interactions central to eCO₂R selectivity. Equivariant graph neural network (GNN) architectures with extended message-passing schemes^[51] partially address this limitation, but a fully satisfactory treatment of long-range electrostatics within the MLFF framework for electrochemical interfaces remains an open methodological challenge, discussed further in Section “Challenges, limitations, and near-future directions of eCO₂R”.

Machine learning regression models

Whereas MLFFs accelerate atomistic simulations by approximating the potential energy surface, a broader class of ML regression models addresses a complementary challenge: rapidly predicting catalytic properties, screening materials, and optimising operating conditions across large design spaces that are inaccessible to systematic DFT computation or experimental screening. Models in use for eCO₂R include decision-tree-based ensemble methods such as gradient boosting regression (GBR)^[53] and random forest regression (RFR)^[54], and graph-based architectures including GNNs^[55,56], all used as surrogate models that map catalyst descriptors, atomic structures, or operating variables to target properties. Decision-tree-based models are especially effective when the input is a tabular descriptor set, as they can rapidly rank catalysts, interpolate structure–property trends, and narrow the search space before expensive first-principles validation [Figure 2]^[57], and are robust to noise and practical for feature importance ranking. Graph-based and transformer architectures become more powerful when the problem requires learning nonlinear relationships directly from atomic environments, adsorption-site geometry, or large combinatorial catalyst spaces^[58]. The principal advantage of these ML models in eCO₂R is therefore not only speed, but scope: they enable a shift from one-candidate-at-a-time DFT computation toward fast surrogate screening, active-learning loops that strategically select the most informative new DFT calculations to maximise model improvement at minimum computational cost, inverse catalyst design, and data-driven optimisation of electrolyser operating conditions^[59].

Figure 2. ML-accelerated prediction of catalytic activity for transition metal Pc DMSCs. (A) Workflow for predicting the limiting potential U_L: DFT- and database-derived features spanning atomic (r₁, r₂), electronic (e_d, N_d, H_f,ox, χ, I_m), and thermodynamic (Q, ΔG_COOH*, ΔG_Max) descriptors are used for model training, model selection, and prediction of target TM₁–TM₂ Pc DMPcs, where TM₁ is the central and TM₂ the side transition-metal atom; (B) Structure of the Pc DMSCs and the transition metals surveyed; (C) Parity plot for the best-performing GBR model: DFT-computed vs. ML-predicted limiting potential (-U_L), with test-set R² = 0.80 and RMSE = 0.12 V. (A-C) Reprinted (adapted) with permission from Ref.^[57]. Copyright © 2021 American Chemical Society. ML: Machine-learning; Pc DMSCs: phthalocyanine dual-metal-site catalysts; DFT: density functional theory; GBR: gradient boosting regression; RMSE: root mean squared error.

Nevertheless, the predictive accuracy of ML regression models depends critically on the quality and coverage of the training data. Most eCO₂R models are trained on DFT adsorption energies computed for clean, idealised surfaces, and their transferability to structurally distinct catalysts, reconstructed surfaces, or co-adsorption environments is often untested. The absence of large, community-standardised benchmark datasets for eCO₂R, analogous to OC20 and OC22 for heterogeneous catalysis^[60,61], represents a significant constraint on model generalisability. Furthermore, descriptor-based models are bounded by the informativeness of the chosen descriptors: if the descriptor set does not encode the effect of surface reconstruction, co-adsorbed spectators, or potential-dependent charge redistribution, the model will produce confident but unreliable predictions for the operating conditions most relevant to eCO₂R. These data and transferability limitations are discussed further as part of the near-term challenges in Section “Challenges, limitations, and near-future directions of eCO₂R”.

Microkinetics, transport, and interfacial thermodynamics

The methods surveyed above describe the catalyst–electrolyte interface at the atomic scale, but several of the phenomena, which are discussed in the catalyst and operando sections that follow emerge only when this atomistic description is coupled to the kinetics of the full reaction network, to mass transport within the diffusion layer, and to the thermodynamics of the electrified interface. Four method tiers bridge this gap and are surveyed briefly here, ahead of the operando discussion that depends on them.

Microkinetic modelling converts the elementary-step free energies and activation barriers obtained from CHE or GC-DFT into steady-state observables, namely surface coverages, turnover frequencies, partial current densities, and Tafel slopes, by solving the coupled rate equations of the reaction network^[62]. For eCO₂R, microkinetic modelling is the tier that links the energetic landscapes of Sections “Conventional fixed-electrons DFT” and “Grand canonical and constant potential DFT” to experimentally measurable polarisation behaviour, and it is essential for rationalising why selectivity between C₁ and C₂₊ products shifts with potential and pH (Section “Electrolyte pH effects”). Its principal limitation is that the rate predictions inherit the accuracy of their energetic inputs, and most existing eCO₂R microkinetic models use CHE-derived thermodynamics, so that the potential dependence of the rates is introduced only approximately, a constraint discussed further as a near-future priority in Section “Bridging the atomistic-to-device scale gap”.

Continuum transport and local-pH modelling describe the diffusion-layer region between the catalyst surface and the bulk electrolyte, solving the coupled transport of ionic species, alongside the homogeneous CO₂/bicarbonate/carbonate equilibria, to predict the local concentrations and the surface pH that the catalyst experiences^[63,64]. Because proton-consuming reduction occurs preferentially at the cathode, these models account for the rise of surface pH by several units above the bulk value under operating current, an effect that is invisible to atomistic models of the surface alone but that strongly governs product selectivity (Section “Electrolyte pH effects”). The principal limitation is that continuum transport is typically decoupled from atomistic detail, and its systematic integration with elementary-step kinetics across the full electrode geometry remains an open multiscale problem (Section “Bridging the atomistic-to-device scale gap”).

EDL theory describes the structure of the interfacial region, the inner and outer Helmholtz planes, the diffuse layer, the PZC, and the electric field across the double layer, using continuum descriptions such as the Gouy–Chapman–Stern and modified Poisson–Boltzmann models, or explicit double-layer DFT^[19]. For eCO₂R, EDL theory provides the conceptual basis for the cation effects examined in Section “Electrolyte cation effects”, where the distinction between inner-Helmholtz-plane coordination of CO₂ and outer-Helmholtz-plane field effects is precisely the distinction between two regions of the double layer. Continuum EDL models capture the mean-field electrostatics efficiently but omit the ion-specific and short-range chemical effects that explicit double-layer DFT and AIMD resolve at considerably higher cost.

Surface Pourbaix (electrochemical phase) diagrams map the thermodynamically stable surface termination and oxidation state as a function of applied potential and pH, by comparing the computed free energies of candidate adsorbed and oxidised configurations under the CHE or GC-DFT formalism^[65-67]. For eCO₂R, this tier predicts which surface state, for example metallic Cu vs. residual oxide or hydroxide, or a given *OH/*H coverage, is thermodynamically favoured under operating conditions, and it therefore provides the thermodynamic baseline against which the potential-driven reconstruction of Section “Applied potential induced reconstruction and property change” and the disorder of amorphous and oxide-derived Cu in Section “Amorphous catalysts” are interpreted. Its principal limitation is that it is a purely thermodynamic construction and cannot by itself capture the kinetically trapped or metastable states that frequently dominate operando behaviour, motivating the dynamic methods introduced above (GCMC, molecular dynamics, and MLFFs).

Multiscale framework and methodological perspective

The computational approaches surveyed in this section constitute a multiscale framework for studying eCO₂R in which each method addresses a distinct but complementary aspect of the problem. Conventional fixed-electron DFT and the CHE model provide thermodynamic free-energy landscapes and activity descriptors for catalyst screening; grand canonical DFT introduces explicit potential control and captures the non-linear electrochemical response of intermediate geometries, binding energies, and surface charge. GCMC extends this picture to fluctuating adsorbate populations and potential-dependent surface compositions through stochastic configurational sampling, while MD-based methods, from AIMD through enhanced-sampling approaches to MLFF-driven dynamics, reveal the real-time interfacial dynamics governing solvent structure, ion specificity, and proton transfer^[68]. MLFFs extend the length and timescales accessible to near-DFT-quality atomistic simulation, and ML regression models enable rapid surrogate screening and data-driven catalyst design across chemical spaces that are inaccessible to systematic first-principles computation. Above the atomistic tiers, the interfacial-thermodynamic and continuum methods of Section “Microkinetics, transport, and interfacial thermodynamics” (surface Pourbaix diagrams, microkinetic modelling, and continuum transport) propagate this atomistic description to the steady-state coverages, turnover frequencies, and local-pH profiles that govern measurable cell behaviour, completing the hierarchy from electronic structure to electrode performance. Rather than competing alternatives, these methods form a hierarchy of approximations whose appropriate use is dictated by the physical question being asked; their capabilities, characteristic outputs, advantages, limitations, and accessible scales are compared in Table 1.

The following sections apply this methodological framework to three catalyst classes, Cu-based materials (Section “Copper-based catalysts”), single-atom and dual-atom catalysts (Section “Atomically dispersed catalysts”), and operando catalyst–electrolyte systems (Section “Operando evolution and co-adsorption”), with the objective of illustrating how the choice of computational approach determines which aspects of catalyst behaviour can be resolved and where the remaining gaps lie.

High-index surfaces

Unlike the ideal Cu(111) and Cu(100) benchmarks, high-index surfaces expose step and kink atoms with lower and more diverse CNs that substantially influence CO₂ activation.

Unlike the ideal Cu(111) and Cu(100) benchmarks, high-index surfaces expose step and kink atoms with lower and more diverse CNs that substantially influence CO₂ activation^[74]. Throughout this review, the coordination number (CN) of a surface atom denotes its first-neighbour coordination, that is, the number of nearest-neighbour metal atoms in the first coordination shell including subsurface neighbours, referenced to CN = 12 in bulk fcc Cu; on this scale, terrace atoms of Cu(111) and Cu(100) have CN = 9 and CN = 8, respectively. The binding strength of CO₂ scales with surface CN: computationally, Cheng et al. showed that the chemisorption energy of CO₂ decreases linearly with increasing CN, and that this trend is amplified under applied cathodic potential, where increased electron density at undercoordinated sites enhances charge transfer to CO₂ and strengthens chemisorption^[73]. Verga et al. demonstrated that 7-coordinated step atoms on Cu(533) exhibit a substantially reduced CO₂ activation barrier relative to the 9-coordinated terrace atoms on the same surface^[75].

Undercoordination also switches the binding mode of CO₂. On Cu(111), CO₂ is physiosorbed through weak van der Waals interactions, while on Cu(100) it is chemisorbed at bridge sites^[76]; the less sterically constrained geometry of stepped atoms instead favours CO₂ adsorption at atop sites, which is associated with stronger binding, greater charge transfer, and more pronounced O–C–O bond bending^[77]. As shown in Figure 3, enhanced CO₂ activation at lower-coordinated step and kink sites reduces the activation barrier of C–C coupling, which is critical for C₂₊ product formation^[78,79]. However, C–C coupling itself proceeds more readily on terrace atoms: Cu(100) terraces exhibit a lower C–C coupling barrier than many high-index step geometries, indicating that step and terrace atoms play complementary roles: step atoms facilitate CO₂ chemisorption and *CO accumulation; terrace atoms govern the kinetics of C₂ product formation. This division of labour is the central mechanistic principle of high-index Cu surfaces, but it also defines their limitation: the periodic step geometry of extended high-index surfaces constrains the accessible coordination environment and cannot reproduce the full structural diversity of finite-size or disordered systems. These considerations motivate the move to nanostructured and amorphous Cu catalysts, discussed in the following subsections.

Figure 3. (A) Profiles of the potential energy surface (ΔG) for the CO–COH coupling on the Cu(711) surface; (B) Schematic representations of the optimised structures; (A and B) Reproduced with permission from Ref.^[78]. Copyright © 2021 American Chemical Society; (C) Coordination number analysis for Cu(911) surface. Reproduced with permission from Ref.^[79]. Copyright © 2022 American Chemical Society.

Nanostructures

Unlike extended surfaces, finite-size nanostructures, including nanoparticles and nanowires, are characterised by a high proportion of atoms at edges and corners with reduced CNs, introducing size- and shape-dependent electronic and geometric effects that substantially alter CO₂ activation and product selectivity relative to planar Cu electrodes. As with high-index surfaces, undercoordinated edge and corner sites show a systematic tendency to strengthen CO₂ binding and activation^[80]. Edge and corner atoms typically exhibit CNs comparable to or lower than those of step atoms on high-index surfaces, and as nanoparticle size decreases, the surface-to-volume ratio increases, exposing a progressively greater proportion of these undercoordinated sites^[81]. However, excessively strong *CO binding at highly undercoordinated sites suppresses C–C coupling by impeding the surface mobility and intermediate proximity required for dimerisation, with the result that C₁ products such as CO, HCOOH, and CH₄ are favoured over C₂₊ species^[81].

At the sub-nanometre scale, the continuous size-effect trend gives way to discrete electronic-parity effects, a consequence of the discrete molecular-orbital structure of sub-nanometre clusters, in which the addition or removal of a single atom switches the ground-state spin multiplicity between open-shell (odd-N, doublet) and closed-shell (even-N, singlet) configurations, fundamentally altering the frontier orbital availability for CO₂ activation. Raju et al. showed that odd-numbered Cu clusters (Cu₃ and Cu₅) exhibited similar behaviour to one another but were clearly distinct from even-numbered clusters (Cu₄ and Cu₆), suggesting that the addition or removal of a single Cu atom can substantially alter activity and selectivity^[82]. Wu et al. similarly showed that breaking the local symmetric structure of Cu₆ clusters promoted the formation of higher-value products^[83]. As cluster size increases beyond the sub-nanometre regime, electronic-parity effects become less pronounced, but the CN distribution and size-dependent effects persist, influencing selectivity. Reske et al. investigated surface CN distribution [Figure 4A], revealing that the ratio of undercoordinated atoms increases dramatically when the diameter is smaller than 4 nm [Figure 4B]^[84]. Mistry et al. compared product selectivity on Cu clusters with different sizes, indicating that ultra-small clusters promote HER while moderate-sized clusters boost CH₄ formation [Figure 4C]^[85].

Figure 4. Size effect exhibited by small Cu nanoclusters. (A) Ball models of spherical Cu NPs with 2.2 and 6.9 nm diameters. Surface atoms are colour-coded according to their first-neighbour CN: CN < 8 (grey), CN = 8 (blue), CN = 9 (red), CN > 9 (green); (B) Population (relative ratio) of surface atoms with a specific CN as a function of particle diameter. (A and B) Reproduced with permission from Ref.^[84], Copyright © 2014 American Chemical Society; (C) Size-dependent selectivity for main products. Reproduced with permission from Ref.^[85], Copyright © 2016 American Chemical Society. NPs: Nanoparticles: CN: coordination number.

Moreover, the geometric shape of nanoclusters independently influences selectivity. Bu et al. reported that icosahedral clusters, exposing predominantly (111)-like facets, showed greater selectivity toward hydrocarbons, while truncated octahedral clusters, presenting both (100) and (111) facets, favoured CO formation^[86]. Nabi et al. found that the most stable pure Cu clusters adopt low-symmetry, amorphous-like geometries, and that the geometric disorder associated with lower local coordination enhances CO₂ activation, particularly for CO production^[87]. The observation that structural disorder at the nanoscale is a feature to be exploited rather than a defect to be minimised motivates the deliberate design of amorphous Cu-based catalysts as a distinct strategy, which is the subject of the following subsection.

Amorphous catalysts

In contrast to ordered crystalline surfaces characterised by well-defined facets, amorphous Cu-based catalysts lack long-range periodicity and instead expose a continuous distribution of local coordination environments, including dangling bonds and undercoordinated sites that are inaccessible on ideal crystalline surfaces. The potential of amorphous Cu for eCO₂R was first demonstrated experimentally by Duan et al., who synthesised amorphous copper nanoparticles and achieved a total Faradaic efficiency (FE) of approximately 59% toward liquid products at -1.4 V vs. RHE, outperforming their crystalline counterparts^[88]. Subsequent work by Xiong et al. reported amorphous CuBi catalysts with Faradaic efficiencies approaching 94.7% for formate^[89], while Chen et al. achieved up to 84% FE toward C₂₊ products using Cu nanoparticles coated with an amorphous shell, with in situ surface-enhanced Raman spectroscopy (SERS) measurements directly linking enhanced *CO surface coverage to the coordination environment diversity of the amorphous phase^[90]. These experimental observations have highlighted the need for computational approaches capable of explicitly capturing structural heterogeneity.

Early theoretical efforts largely focused on partially disordered Cu models. Liu et al. modelled oxide-derived Cu as a partially reduced nanocube containing subsurface oxygen, showing that local Cu distortion strengthens *CO binding and facilitates C–C coupling, though the disorder here is a by-product of oxide reduction rather than true amorphisation^[10]. Song et al. applied AIMD-driven thermal reconstruction of a high-index Cu surface combined with oxygen doping to introduce local coordination disorder and investigate *CO adsorption and dimerisation, though the underlying structure remained anchored to a crystalline parent phase^[91].

To overcome these limitations, Muthuperiyanayagam et al. developed a machine learning force field (ML-FF)- driven melt-and-quench framework to generate amorphous Cu and CuNi surfaces with coordination diversity. Comparison of the crystalline c-CuNi (111) and amorphous a-CuNi (i) surface models [Figure 5A and B] shows that amorphisation produces a structurally disordered surface containing a broader range of local coordination environments. This is reflected in the CN distributions of Cu and Ni atoms [Figure 5C-F], where the amorphous surface exhibits a substantially larger population of undercoordinated surface sites than the crystalline alloy. These amorphous structures were subsequently characterised by DFT and automated adsorption site screening^[92,93]. On amorphous Cu, DFT calculations revealed that undercoordinated sites with CNs equals to 5 and 6 support both CO₂ chemisorption and physisorption, compared to only CO₂ physisorption on the Cu(100) surface^[87], with the chemisorbed state associated with significant charge transfer to CO₂ (-0.48e to -0.98e) and pronounced O–C–O bond bending (123.9-144.8°)^[92]. For C₂ selectivity, amorphous CuNi alloys showed that adjacent undercoordinated Cu-Ni sites (CN = 6) significantly lower the free energy of CO–CHO coupling (ΔG = 0.04 eV) relative to crystalline CuNi (111) (ΔG = 1.57 eV) [Figure 5G and H]^[93]. This stabilisation was attributed to localised charge redistribution and enhanced overlap between Ni d states and O 2p orbitals, enabled by the geometric flexibility of the disordered environment, resulting in a preference for ethanol formation over ethylene^[93].

Figure 5. Atomic models of (A) crystalline c-CuNi (111) and (B) amorphous a-CuNi (i) surfaces colour-coded by CN; Distributions of CN for (C) Cu and (D) Ni atoms on c-CuNi (111); and (E) Cu and (F) Ni atoms on the a-CuNi (i) surface; grey and blue bars denote bulk and surface atoms, respectively; Computed Gibbs free energy diagrams for the electrochemical reduction of CO₂ to C₂H₅OH and C₂H₄ on (G) a-CuNi (i) and (H) c-CuNi (111) surfaces. (A-H) Reproduced from Ref.^[93], under the terms of the Creative Commons CC-BY 4.0 license. CN: Coordination number.

Overall, these studies show that amorphisation introduces chemically distinct coordination and electronic environments absent from crystalline Cu surfaces, providing a mechanistic basis for enhanced CO₂ activation and altered selectivity toward C₁ and C₂ products. Structural disorder is one route to modifying the local electronic environment of Cu active sites; compositional modification through introduction of a second metal offers a complementary and more precisely controllable strategy, which is the focus of the following subsection.

Bimetallic catalysts

Pure Cu catalysts exhibit poor product selectivity and are particularly ineffective for alcohol formation, including CH₃OH and C₂H₅OH. Introducing a second metal dopant generates bimetallic active sites that modify the local geometric and electronic structure, breaking the linear scaling relations between adsorption energies of related intermediates that constrain selectivity on single-metal surfaces and steering reaction pathways towards specific products^[94,95]. Computational studies have established that local coordination geometry governs selectivity more strongly than bulk elemental composition, because asymmetric bimetallic sites interact preferentially with different functional groups of adsorbed intermediates, allowing *CO binding and *CHO or *OH stabilisation to be tuned independently^[94,95]. The identity of the dopant determines the direction of this selectivity shift: noble metal dopants promote C₂₊ product formation, while post-transition metal dopants direct the reaction toward CO or formate.

Noble-metal dopants, exemplified by Cu–Pt, Cu–Ag, and Cu–Au motifs, promote C–C coupling through two complementary mechanisms. First, a *CO supply effect: the dopant acts as a local CO-producing site, increasing *CO coverage on adjacent Cu atoms and the probability of *CO–*CO encounters. Second, a strain-and-ligand effect: the dopant shifts the Cu d-band centre, tuning *CO and *CHO binding energies relative to pure Cu and repositioning the active site closer to the top of the activity volcano. Clark et al. showed that Ag dopants increase *CO surface coverage and simultaneously induce compressive strain at Cu sites, enhancing selectivity toward C₂ alcohols through both effects^[96]. Zhi et al. demonstrated that hydrogen affinity and oxygen affinity of the dopant serve as useful descriptors for product selectivity, with moderate-affinity dopants such as Au and Pd facilitating pathways toward C₂₊ products^[97]. Selectivity between hydrocarbons and alcohols within the C₂₊ manifold is further governed by the oxygen affinity of the dopant toward C₂ intermediates. As shown in Figure 6, dopants with higher oxygen affinity (Rh, Ir, Ru) favour C₂H₄ formation, while those with lower oxygen affinity (Pd, Ag, Au) direct selectivity toward C₂H₅OH. Furthermore, the comparison between RhCu and pure Cu catalysts clearly stresses that Rh doping enhances C₂H₄ selectivity, featuring higher FE and current density, while pure Cu produces more C₂H₅OH throughout the studied potential range^[98].

Figure 6. Electrocatalytic eCO₂R performance for different M-Cu bimetallic catalysts. (A) Oxygen-affinity-dependent selectivity; (B-D) selectivity and activity comparison between Cu and Rh-Cu catalysts. (A-D) Reproduced with permission from Ref.^[98], Copyright © 2024 Springer Nature. eCO₂R: Electrochemical CO₂ reduction; FE: Faradaic efficiency; RHE: reversible hydrogen electrode.

Post-transition metal dopants such as Cu–In, Cu–Sn, and Cu–Zn motifs direct selectivity toward C₁ products through a distinct mechanism. Being more oxyphilic than Cu, these dopants stabilise oxygen-bound intermediates through M–O bonds, facilitating protonation at the carbon atom and thereby promoting formate formation^[99-101]. Zhang et al. showed that trace Sn-doped Cu nanoparticles selectively boost CO formation, while increasing the surface Sn fraction to 50% substantially enhances HCOOH selectivity, a compositional switch that reflects the changing balance between *CO and *OCHO stabilisation as Sn coverage increases^[102]. Wan et al. showed that phase-separated Cu–Zn bimetallic catalysts enhance CO selectivity through a reduced barrier for *COOH formation^[103].

The four Cu-based catalyst strategies surveyed in this section, including high-index surfaces, finite-size nanostructures, amorphous systems, and bimetallic alloys, share a common design logic: each strategy modifies the local coordination environment or electronic structure of Cu active sites to raise steady-state *CO coverage, lower the C–C coupling barrier, or break the scaling relations that constrain selectivity on planar Cu electrodes. Structural heterogeneity, whether introduced through surface steps, finite-size disorder, amorphisation, or a second metal, consistently generates active sites with distinct electronic properties that are not accessible on ideal low-index surfaces and that computational studies based on static CHE-level DFT have only partially characterised. The operando behaviour of these systems, such as surface reconstruction under applied potential, competitive adsorbate populations, and electrolyte-specific interactions, introduces additional complexity that is addressed in Section “Operando evolution and co-adsorption”.

All these designs operate through extended or ensemble active sites in which multiple surface atoms participate in catalysis; a fundamentally different approach concentrates activity at isolated metal centres, maximising atom-utilisation efficiency and offering independent control over the coordination environment and support interactions. This single-atom and dual-atom catalyst design philosophy is the subject of Section “Atomically dispersed catalysts: single-atom and dual-atom sites”.

Single-atom M–N₄ sites on nitrogen-doped graphene

Graphene-supported SACs with N₄-moieties are among the most interpretable model systems in electrocatalytic studies for eCO₂R^[114]. The active site is sufficiently simple for catalytic mechanism investigations and catalytic performance regulation toward specific reactions^[115,116]. Electrically conductive graphene supports that can be readily integrated into electrodes enable efficient electron transfer to atomically dispersed active sites, and support engineering of graphene can tune the local electronic structure and enhance eCO₂R kinetics^[117]. Table 2 lists typical SACs as well as the corresponding main products. Graphene-supported SACs can also exhibit good durability in eCO₂R; for example, isolated Ni single atoms in graphene nanosheets maintained stable CO₂ electrolysis for about 20 h, with post-reaction spectroscopy indicating retention of atomically dispersed active sites^[128]. Unlike conventional Cu-based catalysts that produce diverse products at the cost of poor selectivity, graphene-supported (Fe, Co or Ni) SACs demonstrated exceptional selectivity to C₁ products such as CO, with FE over 90%^{[118,129,130]}. The impressive selectivity toward HCOOH or CH₄ is also reported by altering the metal centre to In or Zn, respectively^[120,121].

The dependence of product selectivity on the metal identity reflects the underlying d-electron configuration of the active centre: computational screening of metalloporphyrin analogues of M–N₄ moieties across 110 systems has shown that d⁶ and d⁷ configurations favour CO₂ chemisorption and charge transfer, while d¹⁰ centres (Cu) are largely inactive toward CO₂ activation, a trend consistent with the selectivity patterns observed in graphene-supported SACs^[131].

Consistent with this trend, Zn and Cu centres in symmetric M–N₄ coordination are weak CO₂ activators and, where active, yield CO as the dominant product. The CH₄-selective entries in Table 2 are therefore not representative of these metals in the N4 motif: methane selectivity is achieved in both cases by departing from it, through the microporous N-doped carbon environment of the Zn catalyst^[120] and through boron substitution to give the Cu–N₂B₂ site^[126]. These results anticipate the coordination-environment strategies examined in Section “Coordination environments of SACs”.

Graphene-supported M–N₄ SACs face an intrinsic limitation imposed by scaling relations: the free energies of *COOH formation (ΔG_*COOH) and *CO adsorption (ΔG_*CO) are linearly correlated across different metal centres because both quantities are governed by the same electronic descriptor, the metal–nitrogen interaction strength, which determines the electron-donating capacity of all eCO₂R intermediates that coordinate through carbon or oxygen to the metal centre^[132]. Therefore, the thermodynamics of CO₂ activation (via *COOH formation) and CO release cannot be independently optimised within the M–N₄ motif by varying the metal alone. A metal centre that stabilises *COOH sufficiently to reduce the thermodynamic cost of the first PCET step tends, along the scaling line, to bind *CO too strongly for facile desorption; conversely, a metal centre that releases CO readily tends to present too high a thermodynamic cost for *COOH formation. This trade-off, which is a statistical correlation rather than a thermodynamic inevitability, explains why most reported graphene-supported M–N₄ SACs are highly selective CO producers, operating near the optimal point of the *COOH/*CO scaling line, but cannot access the further-reduced products CH₄ or CH₃OH, which require sustained *CO coverage and subsequent hydrogenation steps that strong *CO binding precludes^[133].

At the molecular level, NBO analysis of M–POR systems, which are structural analogues of M–N₄ sites, confirms that electron-withdrawing coordination environments increase orbital splitting at the metal centre, suppressing CO₂ activation capacity, while electron-donating ligands enhance it^[131]; this ligand-field sensitivity underpins the rationale for coordination environment engineering as a route to shift the *COOH/*CO scaling line.

Moreover, the dispersed active site cannot stabilise two C₁ intermediates simultaneously, which accounts for the absence of C₂ intermediates^[134]. To break the notorious “scaling relations”, a few promising ideas have been proposed, such as support engineering, synergistic effects in heterogeneous dual atom pairs, asymmetric coordination environments, and ligand co-adsorption^[135].

Alternative supports for SACs beyond conventional graphene

The scaling relation constraints inherent to graphene-supported M–N₄ SACs motivate exploration of alternative supports that can provide less symmetric coordination environments, stronger metal anchoring, and modified charge-transfer characteristics. Graphitic carbon nitride (g-C₃N₄), a two-dimensional polymeric semiconductor with sp²-hybridised C and N atoms arranged in heptazine or triazine units, is the most widely investigated alternative, offering M–N₃ coordination environments that partially decouple CO₂ and CO binding energies and enabling selectivity toward products beyond CO. Although its electrical conductivity is lower compared to graphene, g-C₃N₄-supported SACs offer less symmetric coordination environments, stronger intrinsic single-atom anchoring, and more tuneable charge transfer. Recent studies reported that M-N₃ local environments broadly existed on g-C₃N₄; such asymmetric coordination environments impair the correlation of adsorption energy between CO₂ and CO^[110]; N-rich local environments inherently contribute to stronger anchoring of metal atoms with fewer defects and less local distortion compared with that on N₄-doped graphene, and they provide an effective leverage to tune local electronic structures^[136,137]. More importantly, SACs supported on g-C₃N₄ can steer reactions toward further reduced C₁ products beyond CO. For example, Zhu et al. reported that Cr-g-C₃N₄ boosted CH₄ and CH₃OH formation^[137], while Shui et al. claimed that CH₃OH was the most favourable product on Au- and Ag-g-C₃N₄^[138].

Beyond g-C₃N₄, other two-dimensional materials such as MXenes^[139,140] and MoS₂^[141-144] have also been widely investigated as supports for SACs due to tuneable electronic structures that facilitate HCOOH formation or CO reduction to CH₄ or CH₃OH. In addition to planar 2D supports, carbon nanotubes (CNTs), characterised by adjustable diameters and surface strains, are a typical 1D support for SAC anchoring. Wang et al. reported that sluggish proton-electron transfer was the main bottleneck in eCO₂R, which can be improved by utilising highly conductive CNTs to accelerate reaction kinetics^[145]. Similarly, the curved surface of CNTs also contributes to breaking the “scaling relations”, as previous literature reported that Cr-, Fe-, and Co-N₄ SACs can produce CH₄ or CH₃OH^[146,147], while the asymmetric coordination environments on boron nitride nanotubes (BNNT) open the pathway toward CH₄ and HCOOH^[148-150]. Theoretical calculations indicate that highly curved nanotubes can upshift the d-band centre of the active site^[147]; therefore, retaining adsorbed *CO, while lowly curved nanotubes decrease the interaction with intermediates^[151]. Experimental observations confirm that the FE of CO is negatively correlated with the curvature of nanotubes^[152], which aligns well with theoretical predictions.

Furthermore, it is also noteworthy that phthalocyanine (Pc), featuring a planar N₄-coordinated cavity, is widely used as the support for SACs in electrocatalysis. Recent studies suggest that anchoring the CoPc SAC on highly curved CNTs leads to notable structural deformation and electronic redistribution around the Co active site, which lowers the CO₂ activation barrier^[124]. Intriguingly, other related work claimed that CoPc on CNTs with the appropriate diameter exhibits CH₃OH selectivity, as shown in Figure 8^[153], rather than CO or H₂ as previously reported on Co-N₄-graphene^[118,154]. Since the similarity of the first coordination shell motif between Co-N₄-graphene and CoPc, the distinct selectivity may originate from the discrepancy of global structures. Wang’s group proposed that Co-N₄-graphene and CoPc exhibit different PZC; adding vacancies on Co-N₄-graphene could narrow its PZC gap with CoPc and might facilitate CH₃OH formation on Co-N₄-graphene^[155].

Figure 8. CoPc catalysts supported on a CNT for eCO₂R. (A) Curvature-induced local strain effect; (B) FE of CH₃OH, CO and H₂ in an H-cell with 0.5 M KHCO₃ as electrolyte; (C) FE of CH₃OH, CO and H₂ in a flow cell with 0.1 M KOH + 3 M KCl as electrolyte; (D) Current density in a flow cell. (A-D) reproduced from Ref.^[153] under the terms of the Creative Commons CC-BY 4.0 license. CNT: Carbon nanotube; eCO₂R: electrochemical CO₂ reduction; FE: Faradaic efficiency; SWCNT: single-walled carbon nanotubes; RHE: reversible hydrogen electrode.

The enhanced current density observed in the flow cell configuration [Figure 8D] corroborates the mechanistic picture: curvature-induced electronic redistribution at the Co centre not only steers selectivity toward CH₃OH but simultaneously lowers the kinetic barrier to CO₂ activation, producing measurably higher partial current densities than those recorded on uncurved CoPc or Co–N₄–graphene analogues. Collectively, these findings demonstrate that global structural features, including support curvature and PZC, exert a determining influence on SAC selectivity that is not captured by models restricted to the first coordination shell.

This finding demonstrates that atoms beyond the immediate first coordination exert measurable influence on SAC selectivity, a point that has been underappreciated in studies focused exclusively on the metal centre and its nearest ligands. Recognising that both the support and the global structure modulate catalytic behaviour, the most direct route to active-site engineering remains the deliberate modification of the first-shell coordination environment itself, which is addressed in the following subsection.

Coordination environments of SACs

While the symmetric M–N₄ motif has been the most extensively studied coordination environment in graphene-supported SACs, a range of asymmetric configurations including N₃^[110], N₅^[156], N₃C₁^[157], N₂C₂^[158] and N₁C₃^[125] have been proposed and experimentally validated for eCO₂R, each offering distinct activity–stability trade-offs relative to the N₄ benchmark. Theoretically, Hossain et al. predicted that although the Ni–N₄ motif delivers the highest current at more negative potential, the Ni-N₂C₂ and Ni-N₃C₁ exhibit lower onset potential for eCO₂R^[23]. Han et al. identified the true active sites in Ni–N–C catalysts under operational conditions, establishing that the Ni–N₃C₁ motif is the most robust and active configuration and showing that spatial separation of distinct active sites enhances product selectivity^[159]. Luo et al. systematically substituted N with C in the Ni–N₄ motif and showed that increasing the number of C neighbours strengthens both *COOH and *CO binding, with Ni–N₃C₁ giving the lowest energy cost for the potential-limiting step among the motifs examined^[160].

An important exception to the C-substitution trend is the N₅-coordinated motif, whose special local structure, comprising the symmetric N₄ configuration augmented by one axial N atom bound out-of-plane to the active centre, reduces the *COOH formation barrier and facilitates CO desorption, addressing the principal weaknesses of the N₄ site without incurring the poisoning risk associated with C substitution^[156,161].

M–N₄ coordination confers thermodynamic stability on the active site and suppresses HER through weaker H adsorption, but the correspondingly weak intermediate binding raises concerns regarding insufficient CO₂ activation and high onset potentials for products beyond CO. Substituting N with C atoms improves activity at the cost of reduced active-site stability and an increased risk of intermediate poisoning or elevated HER selectivity, a trade-off that requires careful management in catalyst design.

B-, O-, and S-decorated coordination environments provide further routes to tune catalytic performance. B strengthens intermediate binding and promotes more deeply reduced products: Dai et al. showed that replacing N with B in the Cu–N₄ site strengthens *CO and *CHO binding, promoting CH₄ formation^[126], while Kong et al. demonstrated that B in the first coordination shell promotes CH₄ formation whereas second-shell B has no equivalent effect^[162]. As with C-substituted motifs, B coordination can produce excessively strong intermediate binding, impairing catalytic performance through active-site poisoning. O-coordinated environments substantially weaken intermediate binding, restricting selectivity to CO as the sole eCO₂R product^[127]. As a softer and less electronegative ligand than N, S donates electron density more effectively to the metal centre, partially relieving the excessive metal–N covalency of M–N₄ sites and restoring intermediate binding to a more favourable range; Fe– and Co–S₁N₃ SACs accordingly outperform their N₄ counterparts^[116,163].

Coordination-environment engineering thus provides a powerful and tunable route to control the stability, adsorption energetics, and product selectivity of SACs, and represents one of the most effective strategies for overcoming the limitations of conventional M–N₄ motifs in eCO₂R. A complementary approach, rather than modifying the ligand environment of a single metal atom, is to introduce a second adjacent metal centre, creating a dual-atom site in which cooperative binding between two atoms can stabilise pairs of intermediates simultaneously and enable C–C coupling pathways that are structurally inaccessible to isolated SAC sites; this is the subject of the following subsection.

DACs

DACs consist of dispersed metal atom pairs, either homonuclear (two identical metal atoms) or heteronuclear (two distinct metals). They retain the high atom utilisation of SACs while the second metal centre cooperatively tunes adsorption energetics, charge distribution, and reaction pathways beyond what an isolated atom can achieve^[164]. DACs should not be conflated with the bimetallic Cu systems of Section “Bimetallic catalysts”, which they superficially resemble in containing two metals: the difference is dispersion, not composition. A bimetallic Cu catalyst is an extended alloy whose active ensemble remains many contiguous atoms, whereas a DAC is an isolated, support-anchored pair acting as a discrete two-site unit. Thus, a Cu–Ag pair constitutes a DAC here only when atomically dispersed, not as an extended alloy surface.

Homonuclear DACs serve as a useful conceptual bridge between SACs and heteronuclear DACs: two identical adjacent metal centres introduce dual-site reactivity beyond SAC behaviour while avoiding the additional variable of heterometal compositional asymmetry^[164]. The interaction between the two active sites further refines the local electronic structure and enhances catalytic performance. Zhao et al. demonstrated that the Fe₂N₆ DAC outperformed its FeN₄ SAC counterpart for CO formation by facilitating CO₂ activation and reducing the energy cost of CO desorption at the diatomic site^[165]. Zhang et al. showed that charge transfer within the Pd₂ dimer promotes selective CO formation, with the diatomic character of the active site confirmed by X-ray absorption fine structure (XAFS) analysis^[166].

The neighbouring dual site can also provide sufficient local space and cooperative binding to stabilise two adsorbates simultaneously, which is especially relevant for C–C coupling and C₂₊ product formation. Liu et al. identified FeFe@C₂N as the most active DAC for C₂H₄ production among a computationally screened set, with a CHE-derived limiting potential of U_L = -0.76 V, a thermodynamic estimate of the minimum overpotential that does not account for kinetic barriers along the C–C coupling pathway^[167]. Ma et al. showed that O/OH co-adsorption at the bridge site of Co₂ and Ni₂ dimers destabilises the over-bound bridged-*CO intermediate, a configuration that typically blocks C₂ formation on DACs, thereby opening lower-energy coupling pathways^[168].

Heteronuclear DACs introduce an additional compositional degree of freedom that enables a trade-off between activity, stability, and selectivity inaccessible to homonuclear systems^[169]. The asymmetric coordination of heteronuclear sites allows the two metal centres to favour different intermediate binding modes, C-bound or O-bound, enabling *COOH, *CO, and *CHO binding to be tuned independently rather than being constrained by the scaling relations that limit single-metal and homonuclear surfaces. This asymmetry gives rise to a “division-of-labour” principle in which one metal centre specializes in CO₂ activation while the other stabilizes key intermediates: in the Ni–Fe DAC, for example, the Fe site preferentially activates CO₂ while the Ni site stabilizes *COOH^[170], and one centre may act as a promoter while functioning as a spectator in the direct bond-making steps^[171]. Computationally, heteronuclear DAC performance has been rationalised through descriptors including the *CO binding energy difference between the two metal sites, the degree of inter-metal charge transfer, and the asymmetry of the projected d-band density of states; these are quantities that connect naturally to the ML screening workflows described in Section “Machine learning regression models”. Representative predictions of DACs and their experimental realisation are listed in Tables 3 and 4, respectively. Computational predictions for both homonuclear and heteronuclear DAC systems have been substantiated by experimental realisation. Zhao et al. demonstrated that a precisely synthesised homonuclear Fe₂N₆ DAC delivered higher CO FE than its FeN₄ SAC counterpart, attributing the improvement to facilitated CO₂ activation and reduced CO desorption energy at the diatomic site^[165]. Zhang et al. reported a supported Pd₂ DAC in which charge transfer within the Pd–Pd dimer promoted selective CO formation^[166]. For heteronuclear systems targeting C₂₊ products, Ren et al. demonstrated that isolated Ni–Fe diatomic sites embedded in nitrogen-doped carbon achieved selective CO₂ electroreduction, with operando XAS confirming retention of the diatomic active-site geometry under reaction conditions and synergistic binding between the two metal centres stabilising *COOH^[170]. More recent work has extended DAC synthesis to metal–organic framework supports: Zhao et al. reported that a Cu–Sn hetero-bimetal MOF achieved selective electroreduction of CO₂ to ethanol via an asymmetric *CO–*OCH₂ coupling pathway^[172], while Chen et al. showed that asymmetrically coordinated Cu dual-atom sites on doped carbon directed selectivity towards ethanol with suppressed competing pathways^[173]. Collectively, these experimental results confirm that the cooperative reactivity predicted computationally for DAC systems is chemically realisable, and that the “division-of-labour” principle translates to measurable selectivity advantages under electrochemical conditions.

Despite this progress, DACs present synthesis and characterisation challenges that are qualitatively distinct from those of SACs and that remain incompletely resolved. Achieving precise control over the metal–metal distance and coordination geometry of the diatomic site, parameters that are critical determinants of the cooperative binding effects described above, is considerably more demanding than anchoring isolated single atoms, because diatomic sites must be prevented from either separating into two independent SAC sites or aggregating into metal clusters under the reducing conditions of eCO₂R. Characterizing the true atomic structure of DAC active sites under operando conditions requires a combination of XAFS, scanning transmission electron microscopy (STEM), and ideally in situ SERS, since the diatomic geometry may evolve under applied potential in ways not captured by ex situ structural analysis, a reconstruction pathway that the GCMC and GC-DFT frameworks discussed in Section “Applied potential induced reconstruction and property change” are well positioned to model but that has not yet been systematically applied to DAC systems. Bridging this characterisation gap represents a concrete near-future priority: MLFF-generated structural ensembles for DAC active sites could be used to compute simulated XAFS spectra for direct comparison with operando measurements, providing a quantitative test of whether the computational active-site models correspond to the experimentally observed structures.

Overall, DACs extend the concept of atomically dispersed catalysis by introducing cooperative multi-centre reactivity. Experimental realisations across homonuclear and heteronuclear systems confirm that the adsorption and scaling limitations of conventional SACs can be overcome through rational diatomic site design. The operando structural evolution of these sites under electrochemical conditions, and its influence on product selectivity, is addressed in the context of potential-induced reconstruction and co-adsorption phenomena in Section “Operando evolution and co-adsorption”.

The single-atom, dual-atom, and oxide-derived Cu results surveyed above rest on DFT of doped, heteroatom-coordinated, and partially oxidised surfaces, for which conventional DFT carries limitations more acute than for clean metal slabs. First, exchange–correlation functional error in adsorption energies, of the order of 0.1-0.3 eV on transition-metal and M–N–C centres, is comparable to the free-energy differences that distinguish competing eCO₂R pathways, so that selectivity predictions can be qualitatively sensitive to the functional chosen^[177,178]. Second, on doped and semiconducting supports such as g-C₃N₄ and N-doped carbon, computed adsorption energetics depend on the Fermi-level position and the net charge of the cell^[179]; fixed-electron, fixed-cell DFT represents this dependence only approximately, and doping level, vacancy formation, protonation, and oxidative surface treatment each shift the Fermi level and therefore the computed binding energies of *COOH, *CO, and *CHO^[180,181]. Third, the correlated 3d centres that dominate SAC and DAC chemistry (Fe, Co, Ni) and the partially oxidised Cu of oxide-derived catalysts are susceptible to self-interaction error, which corrupts d-state energetics and spin-state ordering unless a Hubbard-U correction (DFT + U) or a hybrid functional is applied^[182-185]. These limitations are mitigated rather than eliminated: benchmarking adsorption energies against higher-level references, determining U non-empirically (for example by linear response) rather than fitting it, performing hybrid-functional spot checks on spin-state energetics, and adopting grand-canonical or Fermi-level-aware schemes that represent the support charge state explicitly (Section “Grand canonical and constant potential DFT”) each reduce the corresponding error. For the doped semiconductor surfaces on which several of the cited SAC and DAC studies rest, the Fermi-level and surface-treatment dependence of the computed energetics is therefore a determinant of the predicted selectivity rather than a second-order correction, and reported adsorption energies should be interpreted with the functional, U value, and charge treatment in view.

Applied potential induced reconstruction and property change

Simulating catalyst reconstruction and surface roughening under operando conditions remains one of the principal challenges in computational eCO₂R because the active surface can evolve substantially under applied potential and adsorbate coverage^[186]. Using the coupled GCMC–EKMC framework described in Section “GCMC”, Zhang et al. showed that Cu adatoms and small Cu clusters (principally Cu₄) migrate and aggregate at defect sites on the Cu(100) surface under applied cathodic potential^[28], with the resulting Cu₄ clusters positioning adsorbed *CO between the first and second layers of interfacial water, an interfacial geometry that lowers the barrier for *CO activation and C–C coupling^[30]. Cheng et al. demonstrated that simultaneous *CO and *H co-adsorption drives Cu surface roughening by weakening subsurface Cu–Cu bonds through adsorbate-induced charge redistribution, lowering the barrier for Cu atom migration and cluster formation^[73]; the same reconstruction pathway was subsequently identified on Cu(111)^[73]. Together, these results indicate that co-adsorption-driven roughening operates on both low-index Cu facets examined under the studied conditions, evidence consistent with a general roughening mechanism across Cu electrodes under eCO₂R conditions. Beyond structural reconstruction, the applied electrode potential substantially alters reaction energetics, surface coverages, and product selectivity. Alsunni et al. showed that adsorbate geometry is governed by the applied potential, and that the free energy of key intermediates responds non-linearly to applied potential, causing the preferred reaction pathway to switch between moderately and highly reducing conditions; for example, from *CO dimerisation at moderate overpotentials to *CO–*CHO coupling at more negative applied potentials (cf. Section “Grand canonical and constant potential DFT”)^[22]. Using a hybrid neural-network/GCMC approach, Zhu et al. systematically investigated the competition between *CO and *HCOO adsorption on Cu nanoclusters as a function of applied potential [Figure 9]^[187]: *HCOO adsorption, formed via a single proton–electron transfer step from CO₂, dominates the moderately reducing potential window, while *CO adsorption, requiring two proton–electron transfer steps, becomes thermodynamically favoured at more negative potentials. The higher *CO surface coverage at more negative potentials reduces the mean distance between adjacent *CO adsorbates, increasing the density of C–C coupling sites and lowering the activation barrier for *CO dimerisation. Potential-dependent effects are not confined to extended metal surfaces. At molecularly defined and atomically dispersed centres, applied potential modifies the formal oxidation and spin state of the metal, established for CoN₄ sites, where interfacial charge drives an electrostatic spin crossover, and for cobalt porphyrins, where the accessible redox state governs the reduction mechanism^[188,189], redistributing electron density at the metal site and thereby shifting intermediate binding energies and product selectivity.

Figure 9. CO dimerisation on the most stable Cu₈ clusters at different potentials. (A) Distribution of all calculated CO dimerisation barriers on Cu₈(HCOO)₄(CO)₇ at -0.40 V, Cu₈(HCOO)₁(CO)₁₄ at -0.60 and -0.80 V, and Cu₈(CO)₁₅ at -1.00 V. Each line represents an individual CO dimerisation event; (B) Energy profiles of the lowest-barrier CO dimerisation pathways on the four clusters; (C) Transition-state structures correspond to the lowest-barrier CO dimerisation pathways. The bond distances between the two fragments at the transition states are indicated in Å. This figure reprodeced with permission from Ref.^[187], Copyright © 2026 American Chemical Society. RHE: Reversible hydrogen electrode.

Co-adsorption of ligands

Beyond the support and coordination environment effects discussed in the Sections “Copper-based catalysts” and “Atomically dispersed catalysts”, catalytic performance under operando conditions is substantially influenced by co-adsorption of spectator species. These are molecular and ionic entities present at the interface that are not direct reaction intermediates but modify intermediate stabilisation, the local electrostatic field, and the accessibility of reactants and proton donors to active sites.

Because water is the primary proton source in aqueous electrolytes, proton-coupled electron-transfer processes near the cathode generate OH^- locally, leading to high interfacial OH^- concentrations and a significant probability of OH^- co-adsorption at active sites. Sun et al. showed that *OH surface coverage is directly correlated with CH₄/C₂H₄ selectivity on Cu surfaces, with higher *OH coverage favouring C₂H₄ formation^[190]. Further computational studies show that OH^- co-adsorption reduces the *CO–*CO separation distance, extends the *CO surface residence time through electrostatic stabilisation of the adsorbed state, and lowers the C–C coupling barrier^[191-193]. The local electrostatic field generated by co-adsorbed OH^- contributes independently to these effects by polarising adsorbed *CO and *CHO dipoles, a mechanism captured by GC-DFT but absent from CHE-based analyses that treat the interface as field-free.

In addition to inorganic spectator species, organic additives enhance PCET, tune intermediate binding strength, and steer selectivity toward specific products. For Ag-based catalysts, which are already selective for CO production, carefully designed ligand, dihexadecyldimethylammonium bromide (abbreviated 2-C16) co-adsorption increases the FE for CO to nearly 100%^[194]. On Cu-based catalysts, Xu et al. showed that co-adsorbed functional groups exert distinct and group-specific effects on eCO₂R selectivity: NH₃⁺- and SO₃H-decorated Cu surfaces promote CH₄ formation by stabilising *CO, *CHO, and *OCH₃^[195]. Qiu et al. demonstrated the positive role of the NH₃⁺ group for CH₄ formation and the concurrent suppression of HER^[196].

C₂₊ selectivity can also be regulated through tailored ligand modification. Glycine-modified Cu surfaces enhance C₂₊ selectivity via *CO–*CHO coupling, driven by charge transfer between glycine and C₁ intermediates^[197]. A tricomponent copolymer-modified Cu surface achieved 77% selectivity toward C₂₊ products; MD simulations showed that this improvement arises from enhanced reactant access through the porous interfacial structure and from modification of the local electrostatic field^[198].

Co-adsorbed spectator species and organic ligands therefore play an active mechanistic role in tuning intermediate stabilisation, interfacial electrostatics, and product selectivity in eCO₂R, effects that operate in parallel with the influence of bulk electrolyte composition.

Electrolyte cation effects

Alkali metal cations modify the electrochemical double layer, alter the local electric field at the catalyst surface, and, in some cases, directly participate in CO₂ activation, making cation identity and concentration experimentally accessible variables for controlling eCO₂R activity and selectivity.

Two mechanistic interpretations of cation effects have emerged in the computational literature, differing in whether cations act through direct coordination of CO₂ intermediates or through modification of the interfacial electric field. The first proposes that cations directly coordinate CO₂ at the interface and contribute to the first electron transfer step. Qin et al. established that cations switch the operative electron-transfer mechanism: in a cation-free electrolyte, CO₂ reduction proceeds through an outer-sphere electron transfer (OS-ET) at the electrode, whereas alkali metal cations raise the OS-ET barrier and open a competing inner-sphere channel^[19]. Using slow-growth sampling combined with ab initio molecular dynamics (SG-AIMD) on an explicit Au(110)–water interface at U = -1.17 V vs. RHE, Qin et al. showed that, in the absence of interfacial cations, CO₂ activation cannot proceed: the free energy along the adsorption coordinate rises monotonically, and no transition state is located^[19,33]. Introducing two K⁺ ions into the inner-Helmholtz region opens an inner-sphere electron-transfer (IS-ET) channel in which the cations coordinate the bent CO₂^δ- intermediate through short-range K–O interactions, yielding a constant-potential free-energy barrier of ΔG_b = 0.66 eV for the rate-determining CO₂(aq) → *CO₂ step. Under the same conditions, the rate-determining Volmer step of HER carries a markedly higher barrier of ΔG_b = 1.27 eV, so that the K⁺ accumulation at the inner-Helmholtz plane simultaneously promotes CO₂ activation and suppresses water dissociation at the same applied potential, a kinetic blockage effect attributable to disruption of the interfacial hydrogen-bond network at high near-surface cation coverage. Zhang et al. extended this picture through systematic DFT calculations of cation–CO₂ interactions across the alkali series [Figure 10], showing that smaller cations (Li⁺, Na⁺) coordinate only one oxygen atom and activate CO₂ weakly, while larger cations (K⁺, Cs⁺) simultaneously coordinate both oxygen atoms, producing strong activation; larger cations also suppress interfacial proton adsorption, providing a route to inhibit parasitic HER^[199]. Yang and co-workers showed through AIMD simulations that smaller cations are characterised by a rigid hydration shell, causing water molecules to intercalate between the cation and C₁ intermediates and favouring C₁ products; larger cations present a loose hydration shell that creates a hydrophobic microenvironment near *CO intermediates, reducing local proton availability, suppressing *CHO formation, and promoting *CO–*CO dimerisation and C₂₊ product formation^[200].

Figure 10. Local interactions between *CO₂ at the Cu(100) surface and the surrounding environment: (A) Schematic of the bidentate adsorption of *CO₂ at Cu surface in an explicit solvent with C and one of the O atoms coordinating with Cu, and the other O atom forming hydrogen bonds (light blue dashed lines) with surrounding water; (B) Schematic of the side configuration of Li⁺ or Na⁺ coordinated with only one O atom of *CO₂ (magenta dashed line); (C) Schematic of the bridge configuration of K⁺ or Cs⁺ coordinated with two O atoms of *CO₂. (A-C) Reproduced from Ref.^[199] under the terms of the Creative Commons CC-BY 4.0 licence.

A second interpretation proposes that cations influence eCO₂R primarily through modification of the interfacial electric field rather than through direct coordination of intermediates. Ringe et al. showed, using a modified Poisson–Boltzmann continuum model combined with DFT, that cation-induced interfacial electric fields stabilise dipolar intermediates and transition states, including the bent CO₂^δ- species and the *CO–*CO coupling transition state, through field–dipole interactions, accounting quantitatively for the experimentally observed cation-size dependence of CO production on Ag and C₂₊ formation on Cu^[201]. Gu et al. provided further support for this interpretation, showing that physiosorbed hydrated cations near the cathode generate a local electric field that stabilises key intermediates and promotes eCO₂R^[202].

These two interpretations are not mutually exclusive. A reconciling framework distinguishes inner-Helmholtz-plane effects, where cations can directly coordinate CO₂ and stabilise bent intermediates through short-range Coulomb interactions, from outer-Helmholtz-plane effects, where the accumulated cation layer generates the interfacial electric field that stabilises dipolar transition states. Both mechanisms operate simultaneously, with their relative importance depending on cation size, applied potential, and surface charge density. Together, cations reshape the local electric double layer, modify proton availability, modulate HER competition, and in some cases directly stabilise reactive interfacial intermediates. Cation identity and concentration therefore constitute experimentally accessible tools for tuning both activity and selectivity in eCO₂R.

Electrolyte pH effects

Electrolyte pH governs proton availability at the cathode and thereby directly determines the competition between eCO₂R and HER, as well as the distribution of products formed along reduction pathways that diverge at proton-consuming steps.

Because proton-consuming reduction reactions occur preferentially at the cathode, the local pH near the electrode surface rises by several pH units relative to the bulk electrolyte under typical eCO₂R operating conditions. Conventionally, eCO₂R has been performed in near-neutral to alkaline CO₂-saturated bicarbonate electrolytes to suppress HER^[203]. Within the pH range that supports eCO₂R, product selectivity correlates strongly with operational pH. Liu et al. showed, using CHE-based DFT calculations on Cu(211), that higher pH strongly promotes C₂₊ product formation while lower pH facilitates the sequential reduction of C₁ intermediates to C₁ products^[204]. Kastlunger et al. examined the pH dependence of C₂ product selectivity using GC-DFT calculations and showed that the first proton–electron transfer to O atoms of C₂ intermediates produces C₂H₄, which is kinetically preferred owing to the greater accessibility of oxygen lone pairs to proton transfer, while the first reduction at the C atom produces C₂H₅OH, which is thermodynamically favoured but kinetically inhibited under conditions of high proton availability^[205]. Under strongly alkaline conditions, where proton activity is low, thermodynamically favoured C-protonation becomes dominant, a result consistent with experimental observations that C₂H₅OH yield increases substantially at pH > 13.

eCO₂R under strongly acidic conditions has received comparatively less computational attention, owing to reduced catalyst stability and strong HER competition^[206]. Recent efforts have turned to acidic eCO₂R, motivated by the need to avoid the parasitic conversion of CO₂ to carbonate and bicarbonate in alkaline electrolytes, which can consume up to 70% of the CO₂ feed, and by the enhanced ionic conductivity of acidic media. High proton concentrations substantially enhance selectivity toward hydrocarbons such as CH₄ and C₂H₄, but suppressing parasitic HER becomes the central challenge. Three strategies have been proposed: first, the local pH near the cathode rises during eCO₂R as protons are consumed, so that reduction effectively occurs within a neutral or alkaline interfacial microenvironment even when the bulk electrolyte is acidic^[207,208]; second, introducing larger cations such as K⁺ generates a hydrophobic local environment and suppresses proton adsorption at active sites^[209,210]; third, raising intermediate surface coverage and weakening proton binding at active sites further reduces HER competition^[193,211]. The eCO₂R in acidic bulk electrolytes is therefore best understood as reduction occurring within a locally engineered neutral or alkaline interfacial microenvironment, while the acidic bulk conditions serve to minimise carbonate losses and enhance electrolyte conductivity.

The effect of pH in eCO₂R is not merely a bulk electrolyte property but a coupled interfacial phenomenon involving proton availability, the local reaction microenvironment, and competition with HER. The operando complexity documented across he catalyst and operando sections above (dynamic catalyst restructuring, co-adsorption effects, and cation- and pH-dependent interfacial behaviour) highlights a fundamental tension between the idealised models used in most computational studies and the conditions that prevail in real electrochemical systems. Section “Challenges, limitations, and near-future directions of eCO₂R” examines the methodological limitations that underlie this tension and identifies the computational advances needed to resolve it.

Beyond the CHE

As established in Section “Conventional fixed-electrons DFT”, the CHE model, while computationally efficient and valuable for thermodynamic screening, treats each elementary step as a concerted PCET, applies a strictly linear -eU potential correction, and cannot represent the non-linear response of intermediate geometries, binding energies, or surface charge to applied potential^[16]. In eCO₂R, these approximations become particularly consequential because the most mechanistically important steps, CO₂ activation, *CO–*CO coupling, and potential-dependent surface reconstruction, are precisely those where CHE and GC-DFT predictions diverge most substantially, and where the choice of computational model determines not just quantitative accuracy but qualitative mechanistic conclusions. Mechanistically accurate descriptions of eCO₂R therefore require methods that explicitly represent electrode potential, kinetic barriers, and the electrolyte environment.

To address these limitations, grand canonical/constant-potential DFT is a widely used approach to treat electrode potential explicitly. Governed by the constant chemical potential, the number of electrons in a specific system can be varied, which is more realistic than CHE for potential-dependent charge redistribution, double-layer response, and reactions where the applied bias changes the interfacial state^[20].

GC-DFT addresses the potential-representation limitation of CHE by treating the electrode as an open electronic system at fixed electron chemical potential (Section “Grand canonical and constant potential DFT”)^[20]. A productive near-future direction is the systematic integration of constant-potential transition-state methods such as NEB^[212] and dimer^[213] searches within GC-DFT workflows as a standard intermediate tier between CHE-based thermodynamics and full AIMD simulations. Such workflows can locate saddle points and compute activation barriers at a specified electrode potential, enabling assessment of whether thermodynamically favoured pathways identified by CHE remain kinetically accessible under operating conditions. This integration is particularly important for steps where CHE and GC-DFT predictions diverge qualitatively, such as *CO dimerisation at moderately reducing potentials vs. *CO–*CHO coupling at more negative potentials (Section “Grand canonical and constant potential DFT”), and for potential-dependent reconstruction pathways on Cu(100) (Section “Applied potential induced reconstruction and property change”). Coupling constant-potential barrier calculations with microkinetic models would enable genuinely potential-dependent rate predictions, currently unavailable because existing eCO₂R microkinetic models rely on CHE-derived thermodynamics as inputs and would provide a direct computational route to Tafel slopes, partial current densities, and selectivity maps as functions of applied potential, enabling rigorous comparison with experimental polarisation curves.

A parallel limitation is the treatment of solvation. Most eCO₂R DFT studies employ linear continuum solvation models (PCM, VASPsol) that cannot represent the non-linear dielectric response near highly charged interfaces, the specific orientation of interfacial water molecules, or the co-structure of ions and solvent at the inner Helmholtz plane, all of which contribute measurably to adsorption energies and barrier heights in the operando studies reviewed in Section “Operando evolution and co-adsorption”. The development of non-linear implicit solvation models that incorporate field-dependent dielectric saturation, and their systematic validation against AIMD-derived interfacial properties, represents a concrete methodological priority that would substantially improve the accuracy of GC-DFT and CHE calculations without the prohibitive cost of fully explicit solvent treatment.

Bridging the atomistic-to-device scale gap

A further persistent challenge in eCO₂R modelling is the scale mismatch between atomistic simulation and device-level experiment. Most atomistic simulations rely on small, periodic, and idealised interfacial models at the nanometre scale and trajectories spanning picoseconds to nanoseconds, whereas experiments probe catalytic performance across electrode areas of cm², over timescales of hours to days, and under mass-transport conditions (CO₂ diffusion to the electrode surface, local pH gradient formation, and product removal) determined by electrolyte flow rates and cell geometry. This represents a span of roughly six orders of magnitude in both length and time, and descriptors or mechanisms derived from atomistic models under idealised conditions cannot in general be transferred to device-level performance without accounting for the intervening scales of surface coverage, mass transport, and cell geometry.

Recent advances at the atomistic tier are progressively narrowing this gap. Constant-potential GC-DFT provides electrochemically faithful energetics; hybrid explicit/implicit solvation models improve the description of the interfacial dielectric environment; and MLFFs extend accessible system sizes and timescales to scales previously reachable only by classical empirical potentials. A critical prerequisite for integrating these advances into a multiscale workflow is the development of constant-potential MLFFs. These are force fields trained at fixed electrode potential rather than fixed electron count, so that MD trajectories feeding into higher-tier models correctly represent potential-dependent interfacial structure and barrier heights (Section “Machine learning forcefields”).

The most direct near-term route to closing the scale gap is a three-tier multiscale workflow in which MLFF-driven MD provides an atomistic-resolved description of the catalyst–electrolyte interface, microkinetic modelling propagates elementary-step rate constants to steady-state turnover frequencies, and CFD accounts for mass transport, local concentration gradients, and pH profiles across the full electrolyser geometry. Although the individual tiers, microkinetic modelling and continuum transport, are each established (Section “Microkinetics, transport, and interfacial thermodynamics”), their systematic coupling, and particularly the transfer of potential-dependent rate constants from the atomistic tier to the microkinetic tier, remains an open challenge. A further priority is the development of standardised interfaces between tiers, so that atomistic descriptors computed for a given catalyst can be propagated to reactor-level predictions without system-specific reimplementation. A prerequisite for reliable multiscale workflows is also the standardisation of the atomistic inputs that feed into them: the absence of community-curated benchmark datasets for eCO₂R adsorption energies and barriers, across catalyst families, DFT functionals, and solvation treatments, means that microkinetic models built from inputs from different studies are not directly comparable. Development of such datasets, analogous to the OC20 and OC22 benchmarks for heterogeneous catalysis^[60,61], represents a foundational data infrastructure priority for the field.

Treatment of solvent effects

At present, solvent effects in eCO₂R are treated primarily with implicit continuum models such as VASPsol^[214] because they are computationally efficient and compatible with GC-DFT, making them the default approach for electrochemical interface simulations in the eCO₂R literature. Their principal limitation is that the electrolyte is represented as a structureless medium characterised only by its bulk dielectric constant and ionic screening length, leaving atomistic solvent structure, hydrogen-bond fluctuations, proton shuttling, and specific ion or water adsorption without explicit representation. The widely used linear continuum models cannot reproduce the non-linear dielectric saturation and ion-specific response that govern interfacial charging and double-layer capacitance under eCO₂R conditions, and they treat the electrode–solvent boundary as a sharp dielectric discontinuity that excludes charge transfer, exchange interactions, and chemisorption at the first solvent layer^[215-217]. For eCO₂R specifically, this last limitation is consequential for intermediates such as *COOH, whose stabilisation energy depends on orientation-specific hydrogen-bonding interactions with the first water layer that can contribute 0.2-0.4 eV to the adsorption free energy^[215], a contribution that implicit models miss entirely and that propagates into errors in CHE-derived limiting potential estimates. These limitations are most severe precisely where they matter most for eCO₂R: at the highly charged, structured interfacial environment where non-linear dielectric saturation, specific water orientation, and proton shuttling via the Grotthuss mechanism collectively determine adsorption energetics and barrier heights.

Three near-future developments are likely to substantially improve the treatment of solvent effects. First, non-linear and non-local continuum models, including extensions of the VASPsol framework toward field-dependent dielectric response and the implicit electrolyte model of Islam et al., address the failure of standard linear models to capture non-linear dielectric saturation and ion-specific response at charged interfaces^[217]; systematic validation of these models against explicit-solvent AIMD benchmarks for eCO₂R-relevant interfaces, particularly Cu/water and Ag/water systems where MLFF-based trajectories are now accessible^[31,52], is needed to establish their reliability across operating conditions. Second, hybrid explicit-implicit schemes, in which the first one or two solvent layers are treated explicitly at the DFT or MLFF level while the bulk solvent is represented by a continuum, offer a balance between accuracy and computational cost; such schemes have been demonstrated for related electrochemical reactions but remain underexplored for eCO₂R. Third, the interfacial water structures, hydrogen-bond networks, and dielectric response functions sampled in MLFF-driven trajectories (Section “Machine learning forcefields”) can serve as physically motivated reference data to reparametrize next-generation implicit models. Specifically, this data can optimise cavity construction, dielectric permittivity profiles, and ionic screening functions. These quantities are currently fitted to bulk experimental data, leaving them poorly constrained for the highly charged interfacial environments relevant to eCO₂R. Addressing this establishes a parametrisation loop between dynamic and continuum descriptions of the interface, substantially improving GC-DFT accuracy without requiring fully explicit solvent treatment.

Machine learning for eCO₂R

Machine learning represents a promising route for accelerating electrocatalysis simulations, but its robustness, transferability, and generalisation are still limited. For MLFFs, transferability outside the training domain remains challenging, especially when nonlocal interactions, rare reactive configurations, or chemically distinct environments are not adequately represented in the reference set^[45]. More generally, the predictive quality of ML models is still strongly constrained by the quantity, diversity, and consistency of available data. In catalysis, many datasets remain small, biased, or non-standardised, while generating reliable labels from DFT, AIMD, or experiments is itself computationally and experimentally expensive^[218,219]. Another major limitation lies in descriptor engineering: hand-crafted descriptors still exert substantial control over model quality, poor feature choices can constrain performance at an early stage, and a universal descriptor framework for catalytic prediction has not yet been established^[220]. Consequently, model development often remains iterative and system-specific rather than broadly transferable. Recent progress in active learning, learned representations, and data-efficient ML potentials is promising, but these challenges remain active areas of development.

Four near-future directions stand out as having the greatest potential to advance ML-driven eCO₂R research. First, constant-potential MLFFs^[49] represent a critical methodological advance: by training force fields at fixed electrode potential rather than at zero charge, they can capture the potential-dependent reorganisation of the double layer and the associated changes in reaction barriers identified in Sections “Applied potential induced reconstruction and property change”, “Co-adsorption of ligands”, and “Electrolyte cation effects”, without the prohibitive cost of GC-DFT MD. Second, the development of community-curated, standardised benchmark datasets for eCO₂R, analogous to the OC20 and OC22 datasets established for heterogeneous catalysis, would address the data scarcity and inconsistency that currently constrain ML model transferability; such datasets should include not only adsorption energies but also activation barriers, explicit-solvent configurations, and operando structural snapshots. Third, uncertainty quantification (UQ) should become a standard component of ML model outputs in eCO₂R, so that predictions used to prioritise experimental synthesis efforts carry explicit confidence intervals; Bayesian ensemble methods and conformal prediction are practical routes to UQ without prohibitive retraining costs. Fourth, inverse design workflows, in which generative models propose catalyst structures optimised for a target property and active learning loops iteratively refine the training set, have demonstrated early promise for eCO₂R catalyst discovery^[221] and represent the most direct path from data-driven screening to experimental realisation.

Bridging computation and experiment

A cross-cutting challenge that transcends the individual methodological limitations discussed in limitations discussed above in this section is the absence of systematic workflows for integrating computational predictions with operando experimental characterisation data. The classes of catalyst reviewed above are routinely characterised by operando X-ray absorption spectroscopy (XAS), SERS, and in situ transmission electron microscopy (TEM), yet these techniques yield ensemble-averaged or spatially resolved structural information that is rarely used as a direct constraint on computational models. This disconnect means that theoretical predictions of active-site geometries, oxidation states, and intermediate coverages are typically validated only indirectly through comparison with product distributions, rather than through structural observables measured under reaction conditions.

A productive direction is the development of computation-experiment feedback loops in which operando spectroscopic signatures such as XAS edge positions, EXAFS coordination shells, and Raman band positions of adsorbed *CO and *OH are used to discriminate between competing structural models generated by DFT, GCMC, or MLFF simulations. For single-atom and dual-atom catalysts (Section “Atomically dispersed catalysts: single-atom and dual-atom sites”), where the active site geometry is directly accessible by XAS, such feedback is particularly tractable: the approach demonstrated by Martini et al.^[114], combining operando XAS with ML analysis to track active-site evolution in SACs, represents an early example of this strategy and points toward a broader framework applicable across catalyst classes. For amorphous Cu-based catalysts (Section “Amorphous catalysts”), where the absence of long-range order makes conventional crystallographic validation impossible, MLFF-generated structural ensembles could be used to compute simulated XAS or pair distribution function (PDF) spectra for direct comparison with experiment, providing a quantitative measure of model fidelity.

Beyond structural validation, UQ in ML model predictions should be informed to experimental collaborators in a form that directly informs synthesis prioritisation: identifying not the single best-predicted catalyst but the region of compositional and structural space where model uncertainty is lowest and predicted performance is highest. Establishing shared data standards and open workflows for computation-experiment integration represents one of the most impactful near-future investments the eCO₂R community can make.