Realistic Modeling of the Electrocatalytic Process at Complex Solid‐Liquid Interface

Abstract The rational design of electrocatalysis has emerged as one of the most thriving means for mitigating energy and environmental crises. The key to this effort is the understanding of the complex electrochemical interface, wherein the electrode potential as well as various internal factors such as H‐bond network, adsorbate coverage, and dynamic behavior of the interface collectively contribute to the electrocatalytic activity and selectivity. In this context, the authors have reviewed recent theoretical advances, and especially, the contributions to modeling the realistic electrocatalytic processes at complex electrochemical interfaces, and illustrated the challenges and fundamental problems in this field. Specifically, the significance of the inclusion of explicit solvation and electrode potential as well as the strategies toward the design of highly efficient electrocatalysts are discussed. The structure‐activity relationships and their dynamic responses to the environment and catalytic functionality under working conditions are illustrated to be crucial factors for understanding the complexed interface and the electrocatalytic activities. It is hoped that this review can help spark new research passion and ultimately bring a step closer to a realistic and systematic modeling method for electrocatalysis.


Introduction
The urgent need to address the severe environmental crisis and meet the ever-expanding demand for sustainable energy has driven the development of renewable materials and techniques.One of the most promising avenues that can mitigate these environmental and energy dilemmas and rectify imbalances is electrochemistry, [1][2][3] which has been recognized as one of the most thriving means for numerous catalytic processes, to name a few, CO 2 reduction reaction (CO 2 RR), [4][5][6] oxygen reduction reaction (ORR), [7,8] oxygen evolution reaction (OER), [9] DOI: 10.1002/advs.202303677nitrogen reduction reaction (NRR), [10] hydrogen evolution reaction (HER), [11] etc. The success of these electrochemical processes has been evidenced to heavily depend on the screening and rational design of electrocatalysts [12] and effective regulation of the electrochemical systems. [13,14]herefore, advancing our atomistic understanding of the electrochemical environment is of paramount importance.Under realistic conditions, the elementary steps of electrochemical processes usually take place at the solid-liquid interfaces and in particular solid-water interfaces, where the reactants and intermediates not only interact with the catalyst substrate but also dynamically communicate with the solvent molecules.As a result, further developments in these areas require a comprehensive, atomic-level understanding of many phenomena in chemistry, physics, and materials science and provide reasonable guidance for the rational design of electrochemical systems.In electrochemistry, the development of electrode materials remains a significant research focus.Recent decades have witnessed great progress in the development of high-performance electrocatalysts with superior activity, selectivity, and stability, among which, heterogeneous catalysts with atomically dispersed metal sites have gathered considerable attention since their innovative report in 1999. [15]Subsequently, the terminology of single-atom catalysts (SAC) was first proposed in 2011 based on isolated Pt 1 dispersed on FeO x to highly oxidize CO, [16] initiating an intensive exploration of the well-defined active sites of electrocatalysts. [17,18]Downsizing metal sites to the atomic scale endows the catalysts with efficient atom utilization, a unique electronic configuration, and strong metal-substrate interactions. [19][22][23] Other heterogeneous catalysts, such as singlecluster catalysts (SCC), [24,25] double-atom catalysts (DAC), [26,27] and metal-nitrogen catalysts (MNC), [28] have also emerged from the philosophy of constructing well-dispersed metals, forming a prosperous electrochemical community.In addition to developing the now-popularized electrocatalysts, the rational design of active sites, including tuning the coordination environment and electronic structure, is a crucial measure to regulate the overall performance of the catalysts. [29]Besides, metal electrodes have always been flourishing for their incontestable electrochemical performance. [30]he solid-liquid interface is a complex ensemble that requires comprehensive atomic insight to explore the elementary steps of the electrochemical reactions.Enormous efforts, both experimentally and computationally, have been devoted to this endeavor to probe the interfacial domain, [31] which remains a challenging task, even for simple electrochemical systems.The inherent and decorated structures of active sites, along with their corresponding electronic nature are responsible for the electrocatalytic performance. [32]The evolutional interface driven by the applied potential, giving birth to transient metastable states, should be taken into account when modeling the electrochemical process. [33]Aside from the structure-sensitivity character, the detailed properties of interfacial water, such as the physical orientation [34] and hydrogen-bonding network [14] in electric double layers (EDL) have a significant impact on the electrochemical performances of electrode materials.Moreover, the cationdependent interfacial water structure, which yields different electric field strengths, [35][36][37] adsorption rate, [38] stability, [39] and solvation pH, [40] have also been investigated to understand the cation effect on the reaction kinetics.Other complexity arises from the descriptions of the molecular dynamics across the interface, which is highly potential-dependent, and the associated thermodynamic modeling.Despite these challenges, it is necessary to investigate the condensed interface for a comprehensive understanding of electrochemical processes.However, experimental surface characterization methods have limited access to this interface, especially under opening electrochemical reaction conditions.
Although many advanced operando or in situ characterization techniques, such as scanning transmission electron microscopy (STEM), [41] have been employed to probe the electrostatic potential distribution, [42] ion distribution, [43] and water orientation at the interfaces, [34,[44][45][46] the extraordinary complexity of this solidliquid region is such that experimental techniques cannot readily characterize surface active sites and trace the progressive behavior at the atomic level.Additional insight is provided by theoretical studies to unravel the microscopic nature of these interfaces. [47]From a computational modeling perspective, classical molecular dynamics (CMD) simulations have been the most widely employed technique to obtain a microscopic picture of solid-liquid interfaces.Without elaborate thinking, we can cite a bunch of recent examples. [48]Together with coarse continuum theoretical methods, [49] such approaches are capable of connecting microscopic phenomena with meso/macroscopic ones, as they can cover extensive time and spatial scales, typically in nano magnitude.However, they heavily depend on the development of efficient, accurate, and universal empirical force fields and may not provide an exhaustive description of atomistic and electronic structure details, which are indispensable for describing the aspects of interfaces, such as electrochemical kinetics.This is where ab initio molecular dynamics (AIMD) simulations, with their nonempirical treatment of atomic interactions, come into play.AIMD simulations of heterogeneous interfaces represent a promising field, rapidly growing thanks to the advances in the hardware, availability through supercomputing centers, and the development of more efficient algorithms that exploit parallel programming.
In this thematic article, we present recent advances and opening challenges in the computational modeling of the solid-liquid interface, reviewing how various effects are interpreted in the calculations and how they influence electrocatalytic activity and selectivity.First, we concentrate on practical atomistic descriptions of solid-water interfaces, outlining current practices of firstprinciples models for simulating heterogeneous electrochemistry, with an emphasis on our endeavors to realistically consider solvation shells and constant potential.We also present a unifying descriptor for thermodynamically and kinetically translating hydrogen evolution reaction (HER) activities on singleatom catalysts (SACs).In the second session, we discuss critical aspects of the selected well-defined cases.Critical issues here are the structures of the single-atom (SA) sites, including the structure-dependent electronic configuration of the active sites and its developing engineered strategy.Furthermore, we address surface dynamic behavior, interface water, (pseudo-)adsorbates, and metal cations in the interfacial region.How does the applied potential or chemisorbed species drive the reversible or nonreversing restructuring of SAs and metal electrodes?How is the highly dynamic state of the water molecules and metal cations affected by the electrochemical condition?Finally, the upcoming specific challenges in realistic modeling of the solid-water interface will be addressed, highlighting the integration of the theoretical and experimental methods in simulating heterogeneous electrocatalysis to advance their operation in rationally designing a more efficient interface.The field of electrochemical interface science is boundless, and the topics and cases discussed here are inevitably subjective.However, we anticipate that the subjects discussed will constitute important strides toward a more comprehensive understanding of electrochemical interfaces in the future.

Thermodynamic and Kinetic Modeling of the Electrocatalysis
In mechanistic investigations of electrochemistry, density functional theory (DFT) calculations have enabled us to comprehensively understand the reaction processes, which are far from a simple endeavor, but great development has been made in improving models.52] Here, in the following section of this perspective, we specifically focus the subject on the major advances in simulation methods and our recent developments in describing the electrochemical interfaces and electrochemical reactions by utilizing firstprinciple theory.

Computational Hydrogen Electrode Model
For electrocatalysis, one of the most widely employed models is the "computational hydrogen electrode" (CHE), [53] allowing screen over substantial materials from a thermochemical view.In combination with DFT calculations, one can access adsorption energies and the free energy diagram of a specific reaction, taking into consideration of the zero-point energy and entropy corrections. [12]The CHE model regards the electrochemical solvated state of proton and electron as the equilibria of the nonelectrochemical state of hydrogen: 1/2H 2 (g) ⇆ H + + e ˗ , allowing  [57] Copyright 2020, American Chemical Society.c) Selected snapshots from AIMD trajectories show an unconventional dissociative ORR pathway that produces pseudo-adsorbed hydroxide species.d) Probability density distribution of coordination number of O by H (CN O ) in Z-direction, from the three equilibrated AIMD trajectories of *O…OH − , *OH…OH − , and *…OH − , showing the spatial distribution of pseudo-adsorbed OH.Reproduced with permission. [13]Copyright 2022, Springer Nature.e) Statistic numbers of H-bonds between the solvation water and CO 2 reactant along the adsorption reaction coordinate at different potentials on FeNC and snapshots of the solvation environment around CO 2 before and at TS. f) HOMO and g) spin density distribution of the CO 2− anion in a linear or bent configuration.Reproduced with permission. [14]Copyright 2022, American Chemical Society.h) Evolution of the O-H distances during proton transfer and the snapshots of the AIMD trajectories for the *COO s H to *COO a H species on Cu metal surfaces.i) Reaction free energy barrier comparison of CO formation form *COO s H and *COO a H. Reproduced with permission. [58]Copyright 2022, American Chemical Society.j) Charge density difference and k) pDOS and pCOHP comparisons of N 2 adsorption with (lower panel) and without (upper panel) water shell.Reproduced with permission. [59]Copyright 2022, American Chemical Society.
for the determination of the chemical potential for the protonelectron pair.A typical case in point: the DFT-predicted reactivities trend of MNCs toward OER is determined to be Ni > Co > Fe by employing the CHE model to consistently evaluate the Gibbs energies of all species involved in the electrochemical OER network while avoiding the explicit treatment of solvated protons and electrons in an alkaline environment. [12]Effectively, the CHE-based evaluation systems allow qualitatively establishing the correlation between the atomistic structure of the metal centers and the catalytic properties.

Explicit Solvation Shell
[56] Implicit continuum effects are accounted for by utilizing an electric field that modifies the adsorption energies [39] while explicit electrolyte has been manifested to regulate or even subvert the reac-tivity and selectivity of metal electrodes.Recent AIMD and DFT simulations of ORR on MnNC with and without explicit solvation showed that the liquid water environment facilitates charge transfer from substrate to adsorbed *O 2 . [57]Hydrogen bonding stabilization by neighboring water molecules converts the end-on *O 2 intermediate to a side-on configuration (Figure 1a) which is not the preferred geometry otherwise.This transformation with the elongation of the O-O bond distance takes place in a time scale of ≈0.5 ps and stays stable throughout the equilibrated trajectory (Figure 1b).The side-on configuration strongly favors the dissociative pathway of ORR to form *O and the following protonation to form *OH and resulting in a more accurate estimation of the overpotential.By considering a fully explicit solvation scheme, Chen et al. have also discovered an unconventional ORR dissociative mechanism on FeNC electrocatalysts, in which the *OOH intermediate is found to spontaneously dissociate to form a "pseudo-adsorbed" OH − species which spontaneously migrates to the bulk solution phase through the hydrogen bond matrix and being spatially confined at a few water layers away from the catalyst surface (Figure 1c). [13]Such pseudo-adsorbed state in *O…OH − , *OH…OH − , and *…OH − are dynamically confined in one of the shallow, thermal-accessible local minimums on the flat free energy surface dominated by solvation (Figure 1d) so that the dissociated OH − could respond to the redox event at the catalytic center in a coupled manner within a timescale less than 1 ps.Constrained MD investigations reveal that the H-bonds between O in activated CO 2 and the adjacent water molecules are rapidly formed once the TS is reached in the course of CO 2 activation on FeNC (Figure 1e). [14]Such configuration transforms originate from the HOMO-LUMO crossover of the CO 2 moiety, namely, the symmetric nonpolar s-p * HOMO state transforms to the distorted, highly polar p-p * (Figure 1f), which is the original LUMO state, causing the spin density to redistribute to the terminal O atoms (Figure 1g), consequently contribute to a strong H-bond acceptor.Through hydrogen bonding stabilization with an adjacent water environment, CO 2 can be bent and adsorbed on NiNC otherwise will detach from the electrode. [32]In the PCET steps of CO 2 RR, Yan et al. discover the hydrogenation of *CO 2 is more favorable at the O site pointing to solvent (denoted as O s ) on Cu electrocatalysts. [58]However, AIMD simulations suggest an irreversible proton transfer from O s to O a (O a refers to the O of *CO 2 adsorbed on Cu surface (Figure 1h), via the local hydrogen bond network within a short timescale less than 1.5 ps, paving the way for the following dissociation step (Figure 1i).The solvation impact on the electronic configuration of the active center should be reckoned with.By comparing the N 2 activation process with and without explicit consideration of solvent, Qian et al. confirm that the water environment can be responsible for the adsorption of N 2 by facilitating the charge transfer from the metal center of FeNC to N 2 (Figure 1j), furthermore, enhancing the interaction between both  and  pDOS of Fe's d orbitals with N 2 p orbitals (Figure 1k). [59]n extensive AIMD simulation has been conducted to examine the explicit solvation effects on the reaction energetics at a range of charged neutral metal-water interfaces, and compared against the continuum solvent models such as CANDLE. [60]Remarkable deviations of adsorption energies are identified between AIMD predictions and continuum solvation results, particularly for the strongly solvated adsorbates such as *OH and *OOH while no apparent superiority in stabilizing the adsorbates when adopting with and without implicit continuum solvation. [56]The consideration of directional H-bonds and steric water adsorption is evidenced to be essential for an accurate description of solvation at the metal-water interfaces.In conclusion, all these challenges could contribute to the nontrivial character of explicit solvation shells.Notwithstanding, a coin has two sides.The incorporation of explicit solvation inevitably gives us expensive computations.Exploring more advanced algorithms to account for a greater number of explicit solvent environments without compromising ab initio accuracy is the subject of future investigations.

Scaling Relationship and the Volcano Plot
When it comes to thermodynamics, the volcano plots, derived from the scaling relations [61] that linearly correlate the surface binding energies of different adsorbates, including the transition state and reaction intermediate, are one of the most famous reactivity descriptors.The presence of such scaling correla-tion does conveniently compress Gibbs energy diagrams for numerous materials, ranging from transition metals, [61] alloys, [62] metal compounds, [63] and molecular catalysts, [64] and facilitates the simultaneous elucidation of the reactivity trends through the Sabatier-type activity plots, which suggests that the optimal electrocatalyst candidates should possess an optimum binding strength between the adsorbates and surfaces, neither too weak that hinder activation nor too strong that prevent the desorption of the adsorbed species. [65]Recent work by Jiang et al. investigated the relationship between the binding energy of one of the key species *OH (E b, *OH ) and the ΔG of the potential determining steps (PDSs) in CO 2 RR to methane on various Cubased single transition metal atom alloys (TM 1 /Cu (111)) by DFT calculations. [66]The inverted volcano-shaped trend between *CO → *CHO (in blue) and *OH → *H 2 O (in purple) straightforwardly translates the optimal screening of the CO 2 RR electrocatalysts into finding materials with (E b, *OH ) close to the volcano apex (Figure 2e).The V/Cu (111), which is located nearest the volcano apex, exhibits the superior reactivity and selectivity toward CO 2 RR.
Xu et al. have demonstrated that the electrocatalytic activity of MNCs is strongly correlated with the local environment of the metal center, namely its coordination number (CN) and electronegativity, as well as the electronegativity of the nearest neighbor atoms. [67]Specifically, the adsorption energetics of chemically similar adsorbates, such as *OOH, *O, and *OH involved in the ORR, OER, and HER, are found to be highly linearly correlated (Figure 2a).To address this, a universal descriptor  is proposed by taking into account the valence electrons ( d ) of the metal center, combined with the electronegativity of the first coordinated nitrogen and second coordination carbon of the metal center.This descriptor, reputedly, is capable of interpreting the experimentally observed activity trends of the ORR (Figure 2b) and HER (Figure 2c), and outperforming other state-of-the-art descriptors, such as orbital-energy theory [68] and the work function. [69]This  also works for SAC-like metalmacrocycle complexes and effectively reproduces the volcano relationships, which guides further rational optimization of the SACs.However, a subsequent theoretical investigation has revisited the HER activity and evaluated the performance of the  descriptor on a series of MNCs.The results indicate that several factors, such as solvent effects, pH dependence, and even computational parameters of DFT methods, can lead to significant discrepancies in adsorption free energy of ΔG H (Figure 2d). [70]Furthermore, the revised HER volcano curve was obtained by linking microkinetics to the free-energy diagrams, when the applied overpotential and kinetics are taken into account, [71][72][73] highlighting the necessity to consider multiple factors when proposing universal descriptors from DFT screening.
Note, the screening strategy based on the scaling relations trivially assumes that all materials and facets follow the same reaction network, [74] regardless of the reaction site, which inevitably poses a thermodynamic constraint on the studied catalysts and can lead to inappropriate simplification.Our evolving understanding toward scaling relation directs us to break or to circumvent it on account of the oversimplicity of the descriptor when screening or engineering potential electrocatalysts.For example, Chen et al. [75] systematically explored ten heterogeneous doubleatom catalysts for electrocatalytic NRR by DFT calculations and  and c) HER versus the descriptor .Reproduced with permission. [67]Copyright 2018, Springer Nature.d) The comparison results of adsorption-free energies of H between PBE+U and the descriptor  for MNCs.Reproduced under the terms of a CC-BY license. [70]Copyright 2022, The authors, published by ACS.e) Relationship between the reaction Gibbs free energy of the CO, CHO, OH, and H 2 O on TM 1 /Cu (111) as a function of *OH binding energy.Reproduced with permission. [66]Copyright 2022, Tsinghua University Press.f) The linear relationship of the first and the last hydrogenation step of NRR on a series of dual-atom catalysts by DFT calculations.Reproduced with permission. [75]Copyright 2021, American Chemical Society.
found that the positive scaling relationship between the potential determining steps (PDSs) of the first and the last hydrogenation step of NRR, i.e., Δ ads G * N 2 H and Δ ads G * NH 2 compromises the design for efficient NRR electrocatalysts (Figure 2f).As a result, the strategy to weaken the negative impact brought by the scaling relations is adopting optimal electrocatalysts with moderate Δ ads G * N 2 H and Δ ads G * NH 2 (e.g., Cr 2 -N 6 G, Mn 2 -N 6 G, Fe 2 -N 6 G, Co 2 -N 6 G, and Zn 2 -N 6 G), or to break it by the catalyst which favors a low Δ ads G * N 2 H and a high Δ ads G * NH 2 .In addition, synergistic effects in surfaces with various kinds of active sites, such as alloy, allow for circumventing the scaling relationships. [76,77]or instance, PtCu alloy could dramatically accomplish a better catalytic methanol electro-oxidation reaction (MOR) performance compared with the pure Pt electrode. [78]DFT calculations show alloying strategy not only reduces the free energy barriers for binding between *CO and *OH but also contributes to the water dissociation to produce more *OH and the resulting promoted anti-poison reaction, arising from the dual-sites effects.Other breaking strategies by multi-active sites might be achieved using a binary oxide, or promoter-functionalized surface. [79,80]n addition, the volcano plots derived from the Sabatier principle are based on the underlying assumption of the Brønsted-Evans-Polanyi (BEP) relation, which assumes a correlation between the reaction energy and the corresponding barriers.However, it may not hold in cases where the reaction kinetics and thermodynamics are not strongly correlated or decorrelated.For instance, the nonmetallic nature of the SA site on MNCs shouldn't be overlooked in employing the traditional volcano plot to rationally design heterogeneous electrocatalysts for HER. [81]From the Sabatier-principle view, the MNCs with ΔG *H value ≈ 0 show excellent HER rates. [82]The CrNC and FeNC are predicted with the ΔG *H in the range of −0.20 to 0.30 eV. [83]However, the experimentally determined overpotential value of FeNC is ≈0.50 V. [84] The substantial discrepancy appeals for an urgent reevaluation system on SA electrocatalysts.Our very recent theoretical work [85] observes a strongly re-orientated water configuration of *H•••HOH motif (Figure 3b) in which *H•••O-H angle and distance are around to be ≈15°and ≈1.64 Å in AIMD simulations, evidenced by the pronounced peaks in the radial distribution function (RDF) analysis for CrNC, FeNC, and MnNC, located in the *H…H distance range of 1.0 Å-1.8 Å (Figure 3a).The corresponding hydride-like feature in the three MNCs attributes the special neighboring water configuration of *H•••HOH to the considerable accumulated electrons in *H (Figure 3c).Additionally, the abrupt configuration transition of the polarized electrolyte water (Figure 3d) in the course of H adsorption on CrNC by TI simulations echoes the charging behavior by the carbon matrix (Figure 3f) instead of the TM center (Figure 3g) while the gentle charging process on CoNC induces negligible structural disturbance of the surrounding environment.In the case of RhNC, the configurational change of neighboring H 2 O is not observed during the Volmer reaction (Figure 3e).
All the phenomena direct to the charge-dipole interaction (CDI) between *H and neighboring H 2 O, which is quantified by 0 |r i +P j −r j | , leading to the E CDI being briefly translated into the distance (r)-and net charge (q)-dependent kinetics.A universal evaluation strategy toward the Volmer step, derived from the E CDI kinetics in perspectives of coordinated N number, which arises the intrinsic variation of the SA per se identity, [86] integrates not only the thermochemical information of ΔG *H but, more importantly, the kinetics from the net charge  [85] Copyright 2023, American Chemical Society.
of *H views (Figure 3h).This metric re-inspects the kinetically unfavorable region II, where the conventional 'volcano' evaluation fails to describe, and inspires us with the interfacial interaction modulation tactics for the optimal design of the electrocatalysts.This very immediate contribution also emphasizes the critical role of the local solvent environment around the intermediates.'

Potential-Dependent Electrocatalytic Reactivity
The scaling relation derived from the CHE model can facilitate the prediction of optimal electrocatalysts for certain electrocatalytic reactions with fewer independent variables to be considered.Nevertheless, the irrational generalization of protoncoupled electron transfer (PCET) of all the elementary steps in CHE-based models is not straightforwardly applicable to the decoupled proton-electron transfer steps. [87]Sequential protonelectron transfers (SPET) have been prevailingly observed on MNCs [14,59,88] and transition metals [58,89,90] electrocatalysts.The mechanistic scenarios differ strongly across surfaces, applied potentials, and pH. [91,92]To illustrate this point, a complete picture of the activation stages in CO 2 RR on cobalt porphyrin electrocatalyst that allows for the transition between the sequential and concerted proton-electron transfer mechanism is provided by Koper and colleagues. [91]By quantifying the competition between PCET and SPET based on first-principles calculations of acidbase equilibrium constants, they rationalize the experimentally observed pH dependence of CO 2 RR activity (Figure 4a).Recently, a Newns-Andersen model has been reported to determine the rate-determining step (RDS) of the CO 2 RR to CO on TMs, MNCs, and a supported phthalocyanine. [93]Using this methodology, the electron transfer rate can be determined from the width of adsorbate-induced projected density of states (pDOS).Combined with CHE-predicted adsorption energies and pH-dependent activity measurements, ET from TMs to CO 2 to from *CO 2 is determined to be RDS over active potential range, whereas either ET to form *CO 2 or PCET to form *COOH is rate limiting on MNCs (Figure 4b).
In general, such correlations per se do not incorporate reaction kinetics from the effect of applied potential and arbitrarily regard the thermodynamic simulation as potential-dependent kinetics, even without considering barriers which is an essential descriptor in the electrochemical process. [87,89,94,95]In the recently developed promising model of the "constant-potential hybridsolvation dynamic model" (CP-HS-DM) by Liu et al, [96] the competing relation of the formation of H 2 O 2 and H 2 O in the course of ORR on CoNC was investigated.Specifically, they explored the  [91] Copyright 2016, the Royal Society of Chemistry.b) Rate map at −0.8 V SHE and pH = 2 for CO 2 RR to CO across electrodes.Reproduced with permission. [93]Copyright 2021, Springer Nature.c) The competing relation of the formation of H 2 O 2 and H 2 O in the ORR process on CoNC in perspectives of thermodynamic (upper panel) and kinetic (lower panel) modeling.Reproduced with permission. [96]Copyright 2021, American Chemical Society.d) Potential-dependent free energy profiles of CO 2 adsorption at the FeNC-water interface.e) Location of transition states (TS) during adsorption at different potentials.The fitting linear relationship between f) ΔG and ΔG ‡ and g) ΔG/ΔG ‡ versus potential.Reproduced with permission. [14]opyright 2022, American Chemical Society.h) Potential-dependent free energy profiles of CO 2 adsorption (upper) and *COOH formation under basic (lower left panel) and acidic (lower right panel) environment on NiNC.i) LPD-K predicted partial current densities for CO evolution in the basic (left) and acidic (right) phase.Reproduced with permission. [88]Copyright 2022, American Chemical Society.j) Free energy profiles for *N 2 hydrogenation to *NNH on FeNC.Reproduced with permission. [59]Copyright 2022, American Chemical Society.Bader charge evolution in CO 2 chemisorbing process on k) NiNC.Reproduced with permission. [88]Copyright 2022, American Chemical Society.and l) FeNC.Reproduced with permission. [14]Copyright 2022, American Chemical Society.4c).In the kinetic ab initio modeling, another major open challenge of the interfaces is the determination of the driving force, i.e., the electrode potential.Existing methods include the ideal consideration employed in the CHE model.To be more specific, the potential treatment is based on the standard conditions (T = 298 K, pH = 0, p = 1 bar).U = 0 is defined as the electrode potential at which there is an equilibrium between an H + and e ˗ in an aqueous solution and a gaseous hydrogen (1/2H 2(g) ), as stated earlier. [53]Furthermore, the chemical potential of the H + and the e − changes as the [H + ] and the electrode potential varies. [97]Despite the conquer of CHE in electrocatalysis, the explicit influence of varying electrode potentials on, for example, adsorption energies cannot be assessed given that varying electrode potentials lead to changes in the excess charge of the surface and the resulting electric fields.
100] A predetermined number of electrons [98] and counterions of potassium atoms [99] are adopted to establish an approximate double layer.These methods are restrictive in that only an integral and constant number of excess electrons can be induced to the electrode surface, namely, that is under a canonical ensemble (NVT) possessing a varied Fermi level for different adsorbed species in each elementary step.However, since the configurational entropic contribution to the reaction kinetics plays a key role, many excellent studies have evidenced that the NVT MD simulation could yield more realistic solvation-free energies and configurations while not compromising the accuracy too much from the potential variation along the reaction coordinate, based on the premise that multiple constant-charge MDs with wide work function ranges to ensure proper equilibration and extensive sampling. [14,59,88,101]As the thematic issue presented here, AIMD simulations are still on the way to approaching a more realistic scenery of the electrochemical interfaces.
Despite the difficulties in the ab initio treatment of the electrochemical reaction barrier, we do need realistic considerations of kinetics for mechanistic understanding.Elementary-step kinetics of the CO 2 RR mechanism on FeNC and NiNC with a sufficiently thick explicit water layer, [14,88] utilizing the aforementioned methods of surface charge tuning by introducing Na + and K + counterions, are systematically investigated by combining AIMD, constrained MD, and thermodynamic integration (TI).Free energy profiles are constructed at different potentials by TI on the equilibrated constrained AIMD trajectories (Figure 4d).It is observed that both the free energy (ΔG) and free energy barrier (ΔG ‡ ) are linearly dependent on the electrode potentials (Figure 4g) in CO 2 chemisorption courses.The total slope (k = 1.65) of the ΔG-U relationship is not as simple as 1 eV/V on CO 2 adsorption as suggested by the traditional CHE model, originating from the contribution of the potential-dependent solvation effect is actually included during constrained MD simulation.Furthermore, the location of the transition state (TS) along the reaction coordinate shifts closer to the initial state (IS) as the applied potential decreases to a more negative value (Figure 4d).Likewise, this phenomenon is observed in the *CO formation process of *COOH + e − → *CO + OH − on the NiNC surface [102] and the H adsorption step on MoS 2 [103] by CI-NEB methods.Langmuir adsorption model-derived potential-dependent kinetics (LPD-K) [88] is proposed based on the assumption that CO 2 activation should experience an adsorption/desorption equilibrium and the coverage obeys the Langmuir isotherm and can be determined by the applied potential.Using LPD-K formulation, we investigate the potential-dependent free energetics of the CO 2 RR on nitrogen coordinated single-nickel atom catalysts in both acid and basic phase and conclude the RDS under acidic conditions is CO 2 adsorption and under basic conditions, it is the first protonation step (Figure 4h).The calculated partial current densities and onset potentials with 10 mA cm −2 current density by LPD-K (Figure 4i) are comparable with experimental values.Besides the CO 2 RR, the first hydrogenation step of *N 2 to *NNH on FeNC shows potential-dependent characters (Figure 4j). [59]These findings benchmark that with sufficient free energy samplings, the NVT-derived potential-dependent kinetics could quantitatively match and well rationalize the experimental observed electrochemical activity on the MNC electrolyzer.
Using Bader charge analysis, the Fe and Ni centers, however, undergo negligible changes in the MNC motif and remain in their initial charge state.As expected, the charge source for activating CO 2 is mainly contributed by the charged Ndoped graphene substrate (as an electron reservoir) (Figures 4k,l).As aforementioned, each constrained AIMD simulation is performed within the canonical ensemble with constant charge, which consequently comes with a shift in the work function along the reaction coordinate and the resultant varied electrode potential.Even with the adopted potential value (U r ) of the studied elementary step being the average one of U IS and U FS (i.e., U r = (U IS + U FS )/2), strictly speaking, the methodology is not constant.Obtaining an accurate free energy landscape at a constant potential, i.e., doing sufficient sampling for both the configurational entropy and electronic contributions within the grand canonical ensemble (of electrons), has been remaining a challenging task.

Constant Potential Models
Theoretically, one of the most widely implemented firstprinciples DFT methods to achieve the constant applied bias potential is fixing the work function or Fermi energy level of the electrode materials by adapting the number of electrons, n, to the applied constant potential, and then being referenced to an absolute potential of an electrode, such as SHE. [104]The implementation of implicit solvation models enables continuous variation in the number of electrons of an electrolyte. [105,106]As a consequence, the grand free energy of an electrochemical system at a fixed bias is obtained by variationally optimizing n, and the generated charge will be balanced by ionic screening in the electrolyte.For example, Xiao et al. used the implicit CANDLE [60] solvation model implemented in JDFTx [106] to investigate the atomistic mechanisms underlying electrochemical CORR on Cu (111)-water interface, which including both the constant potential and solvation effect. [105]These simulations successfully elucidated the competition among various pathways of CORR on Cu at different pH and accurately predicted the onset potentials from QM calculations, validating the particular combined consideration of constant potential and implicit electrolyte effect.It's worth noting that the introduction of n as an additional variable causes numerical instability for the self-consistent field (SCF) method when solving the Kohn−Sham equations.The very recent publication by Xia et al. has formulated an efficient and robust fully converged constant-potential (FCP) algorithm based on Newton's method and a polynomial fitting to accelerate the convergences, outperforming other algorithms, such as CG and BFGS. [107]Benchmarked constant potential MD simulations of the first electrochemical hydrogenation step of CO catalyzed by the FeNC at U = −1.5 V revealed that the hydrogenation step is tend to occur via the CHO pathway, in agreement with previous experimental reports, exhibiting the accuracy of the grand canonical ensemble (μVT) modeling of electrochemical interfaces for the quantitative estimation of macroscopic electrochemical observation.
Although fractional charges, being introduced to achieve continuous work function and the corresponding controllable potential, [100] offers the chance of electrochemical modeling at a constant potential since the Fermi level can be matched to the desired bias iteratively.The implicit solvation also brings about some extra errors and artifacts, for example, the induced balanced charge is unphysically dispersed in the electrolyte and renders the incapability of capturing the impact of double-layer charging on reaction energetics.The use of work function as a descriptor of the driving force leads to cell size-dependent functions of reaction energetics.In addition, the ignorance of explicit solvation shells, as aforementioned, will necessarily hinder our mechanistic knowledge as well.
An alternative to realizing the constant potential simulation is transforming the energetics based on the charged neutral model (CNM) to that on the constant potential model (CPM).Cell [108] and charge [109] extrapolation schemes based on fully explicit treatments of the electrolyte to deduce CPM, which are attainable for the evaluation of kinetic landscapes, both of which accommodate the effective surface-charge density to appropriately depict the electrostatic effects of the double layer on reaction energetics [110] on account of the description of the variations of the interfacial field local to the active region.However, recent studies [52] have shed light on the complexity of the proposed cell-extrapolation scheme, such as the inevitably prohibitive cost to apply to more complex electrochemical processes and larger-size supercells.Due to the assumption that the "chemical" and electrostatic contributions to the energetics are separable (i.e., E = E chem + E el ), the Figure 5. a) GCP-K predicted Tafel plots for HER on MoS 2 electrode in both acidic and basic conditions.Reproduced with permission. [103]Copyright 2018, American Chemical Society.b) The experimental and theoretical quantum capacitance of pristine and N-doped graphene as a function of electrode potential.Reproduced with permission. [115]Copyright 2012, the Royal Society of Chemistry.
charge-extrapolation which defines that the interface functions as a simple capacitor, contradicts the behavior of the physical systems, especially in the case of 2-dimensional materials.
Grand canonical potential kinetics (GCP-K) [102,103] are derived from minimizing the fixed-charge free energy, F(n), to grand canonical, G (n, U) using Legendre transformation, allowing both the geometry of the transition states and the charge transfer from heterogeneous electrodes to adsorbates to change continuously along the reaction coordinate as the potential is changed.Using this GCP-K predicted free energy, the hydrogen evolution reaction (HER) at a sulfur vacancy on the basal plane of MoS 2 in both acidic and basic conditions is predicted. [103]The higher activity of HER in basic originates from the TS of the Volmer step locates closer to FS compared with that in acid, leading to a more favorable Tafel slope of 60 mV/dec (Figure 5a).
We notice that a quadradic term of F (n) is included in GCP-K formulations, specifically, in which GCP (U) depends quadratically on the number of electrons, n, and the applied potential U, allowing a continuous description of the evolution of transition states.This quadratic grand canonical potential GCP (U) accounts for the change in capacitance as the potential changes.However, we cast doubt on the proposition because charge variation on 2D materials (such as the adopted representative cases of N-doped graphene [102] and MoS 2 103 in GCP-K) can have a much stronger impact on the electrochemical reaction than the charge on 3D metals, [111] without mentioning the widely employed fixed quantum capacitance value (C d ) of 21 μF cm −2 of single-layer pristine graphene deviate from experimental values. [112]Way back in the decade, experimental [113][114][115] and theoretical [116] scientists have reported that the quantum capacitance of C d is linearly dependent on the absolute potential for both pristine and N-doped graphene electrodes (Figure 5b).The unconformity between the realistic potential-dependent quantum capacitance and the adoptive constant one in the electrochemical elementary process will inevitably arise nonnegligible variations in accessing the reaction behaviors and stand in the way of the realistic interface.On the other hand, from the implicit solvation perspective, the adoption of continuum approaches in conjunction with explicit solvation shells in CP-HS-DM [96,117] offers a more realistic description of hydrogen bonding and solvation.
To the best of our knowledge, there was no theoretical continuous potential described elementary reaction thermodynamics and kinetics while the experimental investigations concentrate on the overall kinetic performance under continuous potential conditions.The linearity between the capacitance and potential, instead of the traditional consideration that the surface charge density linearly correlates with potential, calls for a pressing re-examination of the potential-dependence property.On the other hand, emerging hybrid explicit-implicit approaches have the potential to advance the field of constant potential modeling, since they account for the solvation with cheap computational cost while retaining explicit electrolyte, representing a cuttingedge method. [110]However, developing effective implicit models that accurately describe the smooth transition of the permittivity between implicit and explicit solvation is another consequent issue.

Structure and Dynamics of Electrocatalysis at Complex Solid-Liquid Interfaces
Distinct from the individual bulk media, solid-liquid interfaces, and particularly solid-water interfaces play significant roles in various territories of modern chemistry, such as heterogeneous electrocatalysis, and new materials design.To gain a deeper understanding and promote applications in a wider domain, a detailed comprehension of the interfaces including physical configuration, chemical and electronic structure and response to the environment at the atomistic level is a prerequisite to grasping the electrochemical process and its performance.In this section, we will specifically focus on the molecular representation of the electrode and water, and accordingly, examine the structure-activity relationship.An elaborate microscopic description of solid-liquid interfaces can be obtained from a theoretical point of view by AIMD simulations, which will be the primary approach used in this section.At this atomistic and molecular scale, the key questions to be addressed are, i.The molecular engineering at the interface from the perspective of active sites;  111), Pt doped @WC x , and Pt ads @WC x .d) Gibbs free energy of *H adsorption for different catalysts.e) The Bader charge analysis of the atomic Pt site.f) PDOS of Pt 5d orbitals and W 5d orbitals in the Pt doped @WC x and Pt ads @WC x .Reproduced with permission. [120]Copyright 2022, Wiley-VCH GmbH.g) Free energy diagram for ORR reaction pathways on Cu 1+ -SA and Cu 2+ -SA active sites in alkaline media.(pH = 13) h) PDOS overlap for Cu(I) and Cu(II) with *O atoms.Reproduced with permission. [121]Copyright 2020, American Chemical Society.i) Linear correlation of the CO 2 activation barrier, ΔG ‡ , versus U RHE on NiNx electrodes.j) Bader charge of Ni center for NiNx under PZC.k) PDOS of 3d x2-y2 orbitals for Ni center in all sites and C 2* orbitals of CO 2 .Reproduced with permission. [88]Copyright 2022, American Chemical Society.
ii.How the chemical and physical configuration affects the electronic nature of the sites and the related electrochemical performance; iii.The dynamic evolution of the potential-and pH-dependent adsorbate coverage and the corresponding induced restructuring of electrodes and adaptive behavior of water molecules; what's more, we emphasize the dynamic behavior of reaction species, i.e., the dynamically confined 'pseudoadsorption' catalytic state; iv.Considering the complexity of water solvation in realistic electrochemical systems, where ions are present at various concentrations, the effect of cations and interfacial pH must be incorporated since they are key elements that influence the relative organization of both the solid and the liquid as they come into contact continuously.
The structure-performance relations of these active sites, their dynamic responses to the environment, and catalytic functionality under working conditions are crucial factors for understanding the interfaces. [118]

Engineering Electrocatalysts Via Structure Regulation
Heterogeneous electrocatalysts are intricate structures that are typically composed of several phases, generally, including an active catalytic phase and a supporting one with a conducting matrix onto which the nanoparticles of the active phase, such as transition metal atoms are embedded which not only stabilize but also collaborates with the active sites. [119]Platinum-based electrocatalysts occupy a pivotal position in the diverse catalytic process, but their scarcity and high cost have hindered large-scale application.In our recent work, the crystalline lattice-confined atomic Pt in tungsten carbide (Pt doped @WC X ) exhibited competitive HER performance with Pt@C. [120]Notably, DFT calculations revealed that the existence of the support phase of WC decreases the water dissociation energy barriers (Figure 6b), leading to 40 times greater mass activity than that of the commercial Pt@C in alkaline phases.Charge density differences disclose that Pt doped @WC X exhibits near-zero valence states (−0.14 |e|) and matched electronic structures as metallic Pt (Figure 6c), delivering similar catalytic behaviors when adsorbing H in the acidic phase (Figure 6d), which are much superior to the adsorbed and NiFe-DASC (lower).b) Calculated free energy diagrams for CO 2 RR and OER processes.Reproduced with permission. [29]Copyright 2021, Springer Nature.c) Proposed reaction scheme of the single-site (upper) and dual-site (lower) mechanisms towards OER.d) Calculated free energy diagrams at 1.23 V for OER over FeNC, CoNC, and NiNC.Reproduced with permission. [12]Copyright 2018, Springer Nature.e) Calculated free energy diagram of NiPc-CN-H 2 , NiPc-H 2, and NiPc-OMe-H 2 electrocatalysing CO 2 RR at −0.11 V SHE .Reproduced with permission. [124]Copyright 2017, Springer Nature.f) Electron density difference plot showing charge transfer from graphene substrate to NiPc molecules.Reproduced with permission. [126]Copyright 2021, American Chemical Society.
Pt sites on WC (Pt ads @WC X ), negatively charging by −0.65 |e| (Figure 6e).The different coordination pattern between the active site, Pt, and support matrix, WC, arises a huge deviation of the electronic configuration (Figure 6f), thus of the electrocatalytic activity.The novel electronic character of the SA site and the heterogeneity of its local environment cast doubt on mechanistic studies.The support of WC can also be a direct participant in the reaction.Despite both the water dissociation barriers on Pt ads @WC X (0.35 eV) and Pt doped @WC X (0.15 eV) are very low compared with that of metallic Pt (111) (0.86 eV) in the alkaline phase (Figure 6b), their dissociation properties are not identical, especially, regarding the adsorption site of produced *H.To be specific, on Pt doped @WC X , *H species adsorb directly at the Pt site, while on Pt ads @WC X , they adsorb at the W site, resulting in an additional 0.47 eV free energy compensation for *H migration from the W site to the Pt site during the subsequent H 2 generation.
The character of structure-sensitive performance is embodied in MNCs as well.Engineering the solid-liquid interface by mixing the Cu salt precursor with varied urea content, a tailored coordination environment of the Cu (I) site with corresponding boosted ORR activity has been reported by Sun and coworkers. [121]Theoretical calculations confirm the lower CNs with neighboring N of Cu feature a lower oxidation state, which facilitates the RDS of *O 2 protonation to *OOH proceeding.The CHE-predicted ΔG is lower on Cu (I) site by 0.25 eV compared to Cu (II) site (Figure 6g).The reactivity difference could be attributed to a less occupied O 2p-d antibonding state in Cu (I), as revealed by the pDOS in Figure 6h.In the electrocatalytic conversion of CO 2 to CO by NiN X ( X = 1-4) moiety, our erewhile theoretical investigations by DFT-based AIMD simulations confirm coordinately unsaturated nickel-nitrogen sites on graphene (NiN 1 ) achieve higher current density than others. [88]Particularly, in the lower electrode potential regime of U RHE > U 1 , NiN 1 , and NiN 2 are more active than NiN 3 and NiN 4 towards CO 2 RR (Figure 6i), stemming from the lowest oxidation state of Ni in NiN 1 (Figure 6j) and its occupation of the 3d x2-y2 orbital (HOMO) of the active center, Ni, which is more inclined to overlap with the C 2* (LUMO) orbital of CO 2 (Figure 6k).Such catalyzing trends are ascribed to the charge capacity of Ni when the NiN 1 structure has the lowest barriers for CO 2 RR and a high barrier for the competing reaction HER. [32]In short, the modification of the metal SA with a well-controlled coordination environment and the resulting modified low oxidation state can significantly contribute to the performance of SA sites.
NiN 4 moiety can serve as a structural promoter (Promoters are typically hetero-atoms or hetero-moieties that spread over the active surface to enhance catalytic activity or selectivity) of CoN 4 and FeN 4 sites in diatomic sites (DACs) to promote electrocatalytic performance. [29,122]Combined DFT and AIMD calculations reveal that the adjacent NiN 4 site can effectively adjust the electronic localization of the proximity CoN 4 site, promoting the *OH desorption and *H adsorption on the CoN 4 site. [120]Simultaneously, the orbital coupling between Fe and Ni lowers the molecular orbital energy levels of Fe and adsorbates (Figure 7a) and weakens binding strength to the reaction intermediates, thus boosting CO 2 RR and OER performance (Figure 7b). [29]It is worth noting that the active sites are not limited to the metal center for MNC electrocatalysts.Huang, Duan, and coworkers analyzed the relationship between the atomistic structure of the metal centers of Fe, Co, and Ni moieties and OER activities. [12]It is theoretically predicted NiN 4 exhibits the highest performance as revealed in Figure 7d, because of the collaborative participation of C atoms (a dual-site mechanism, Figure 7c) in the OER process.
However, these 2D materials suffer from the instability of the sintering of the nanoparticles of the active component with a low surface area. [123]Additionally, the uncontrollable doping environment leads to difficulties in defining the local electronic structure of the SA sites.To address these issues, researchers have explored and designed a new class of molecular electrocatalysts.In 2017, ZHANG et al. reported the noncovalent immobilization of molecularly engineered cobalt phthalocyanine on carbon nanotubes (CoPc@CNTs), achieving superior faradaic efficiency (FE) for CO 2 RR to CO. [124] The structural uniformity and tunability, as well as the - interaction that preserves the molecular integrity and fine-tunes the electronic structure of the active sites.Subsequently, a series of further engineered NiPc@CNTs achieved unit conversion to CO with high current density, with tetramethoxy substitution as an electron-donating group (EDG) improving the stability of NiPc, and cyano group (-CN) as an electron-withdrawing group (EWG) destabilizing the catalysts in long-term electrocatalysis. [125]Atomic charge analysis revealed that the introduced -OME to the NiPc enriched the electron density of the NiN4 moiety, resulting in strengthened COOH binding under the CHE scheme, while weakened CO adsorption prevented the poisoning of the active site (Figure 7e).Furthermore, the Ni-N bond was strengthened by the more electron-rich metal center, which protected the pristine NiN4 moiety from dissolution and the resultant deactivation. [126]In our recent theoretical study, the interaction effect and charge transfer kinetics between NiPc molecules and the nanocarbon support were investigated using the electron density difference map. [126]The enhanced polarization of NiPc-CN and the resultant stronger interaction with the graphene substrate (Figure 7f, upper panel) facilitates the charge transfer across the interface, which is significant in electrocatalysis but is often neglected in computational modeling.These findings open up an additional dimension in molecular design.

Dynamics of the Solid-Liquid Interfaces
Hitherto, a plethora of experimental and theoretical reports have been focused on the various electrocatalysts model to investigate their mechanistic behavior and reactivity trends.However, solely understanding the interface in its resting state may not always be sufficient to rationalize the electrocatalytic performance.Our findings suggest that Our obtained knowledge and findings suggest that our scope should not be confined to stationary form. [17,18]During catalytic conditions, however, the surface of heterogeneous electrodes, together with the water molecules are subject to surrounding mediums, i.e., the electrochemical potential, pH, temperature, and the coverage of adsorbates, which could render the catalytic electrode dynamic, make the surface reconstruction, [127,128] and induce transient metastable surface species. [75,129]Capturing and harnessing the dynamic distribution of the direct participating species at different physical scales in the electrochemical interface is of importance when modeling the active interface.This section specifically highlights this particular research field: the dynamic behavior of the interface, which incorporates the leaching behavior of the SA sites and surface restructuring of the metal electrode; the adaptive coordination and orientation of the water molecules; and the dynamically confined 'pseudo-adsorption' catalytic process.Advancements in theoretical understanding make it possible to harness the dynamic interface to a greater extent and in a more systematic way.

The Leaching of the SA Sites
The electrocatalytic performances of electrode materials are inherently determined by the number and structure of the active sites, which are critical and yet less controllable due to the complexity involved in the electrocatalytic process.132] For example, the genuine site structure for ORR in the case of CuNC electrocatalyst has ever been claimed to be CuN 2 , [133,134] CuN 3 , [135] and CuN 4 , [136] without a consensus arise not only from the different preparation methods, which may result in diverse active structures but also from the challenge that the transient lifetime of SA sites and the limited spatial resolution makes it difficult to probe directly in atomistic dimension by experimental methods under operating conditions.In fact, the dynamic evolution of the active site is ubiquitous in catalysis. [137,138]Wang et al. have studied the reaction mechanism of CO oxidation on ceriasupported Au catalysts static using DFT calculations and AIMD simulations. [18]An initiative dynamic single-atom catalytic mechanism was identified, in which CO adsorption engenders an isolated cationic Au + ˗CO unit's formation first, facilitating the subsequent CeO 2 reduction and CO oxidation.After the catalytic cycle, the Au single atom is reintegrated into the Au nanoparticle, indicating that the actual catalytic single-atom species only exists as a transient under operating conditions which may not necessarily be easy to capture ex situ.Besides, the adsorbate-induced novel dynamism is highly relevant to both the temperature and oxygen partial pressure. [17][141] ZHAO et al. utilized CP-HS-DM [96] to investigate the dynamic progression of CuNC under reducing potentials. [117]The highly thermodynamically favored H adsorption on N sites in the large potential regime of U RHE < −1.0 V, as shown in Figure 8b, is the direct driving force of the notorious restructuring phenomenon of CuNC.After adsorbing one or two H + , two transient states of Cu, with the incomplete leaching being tethered to one N atom (Figure 8c) and the complete leaching dissolved in water (Figure 8d), have been identified at a dynamic electrochemical interface by MD simulations within a timescale <0.5 ps.The reversible transformation between Cu SA and Cu 3/4 clusters is further declared by the MD observations, in which the two leached copper atoms approach one another and concomitant by the elongation of the Cu-N bond and in the thermodynamically favorable dispersion process, the corresponding Cu-Cu bond elongation is identified.Be that as it may, subject to the expensive computational cost, the following aggregation of the simple Cu 4 clusters and the reverse fragmentation process is not completely depicted.This work highlights that the adsorption of H, depending on the potentials, is a crucial driving force for Cu SA to Cu clusters, while the reverse process is dominated by Cu oxidated by OH radicals.As far as we concerned, this work is far from satisfactory on account of the discrete potential simulation and arbitrary introduction of the He atom and OH radicals.In addition, they unreasonably concluded the processes of SA aggregation, CO 2 RR, and decomposition of Cu 4 clusters are sequential (Figure 8a).The considered reducing potential for the CuNC reconstruction, −1.2 V, is in the CO 2 RR active region [142,143] so that the evolving Cu and the resulting undercoordinated configuration are never static but accompanied by electrocatalyzing the CO 2 RR.As many literatures pointed out  and d) the incomplete one.Reproduced with permission. [117]Copyright 2022, American Chemical Society.e) Calculated current densities for the Pt (111), Cu (111), and the undercoordinated Cu adatom active site motif generated at the restructuring Cu electrode.Reproduced with permission. [152]Copyright 2020, Springer Nature.f) Calculated surface Pourbaix diagram showing potential and pH dependence of the Cu (111) surface state.g) Reactive evolution of the A1-15 structure.Reproduced with permission. [127]Copyright 2023, Wiley-VCH GmbH.h) Calculated surface Pourbaix diagram showing potential and pH dependence of the Cu (100) surface state.i) Total grand canonical free energy as a function of H coverage at different electrode potentials.j) Color-filled contour map of Cu (100) with 1 H showing a weakening of the Cu-Cu bond between the top layer and the sublayer.Reproduced with permission. [128]Copyright 2022, American Chemical Society.
that the transient low-coordinated Cu SA sites have a much better performance toward CO 2 RR. [144,145]As a result, we wonder about the genuine active site for the presentative high electrocatalytic CO 2 RR on CuNC under operating conditions.In addition to theoretical simulations, numerous experimental explorations identified the reduction of Cu 2+ to Cu + and Cu 0 and the subsequent aggregation of Cu 0 single atoms occurs concurrently with the enhancement of the NH 3 production rate and ORR, both of them are driven by the applied potential switching from 0.00 to −1.00 V versus RHE. [10,33]By combining operando XANES spectroscopy and DFT calculations, the potential driven dynamic evolution of CuNC from Cu 2+ -N 4 to Cu + -N 3 and further to HO-Cu + -N 2 has been unambiguously identified. [33]his transformation of the Cu SA site is ascribed to the fact that the more negative potential drive more electron aggregate on the catalyst surface, promoting the adsorption of H + , which is a prerequisite driving force for the leaching of Cu single atoms from the catalyst surface and transforming into Cu small agglomerates.From the micro-electronic perspective, the fully occupied 3d orbitals of Cu atom [146] and the consequential weak Cu-N bond, induced Cu and N sites tended to react with H + and resulting in the breakage of the Cu-N bond.This is the underlying reason that the leading actor in the leaching phenomenon is notorious copper.However, NiNC, [147] CoNC, [13] and FeNC [13,148] echoed.

Surface Restructuring
The underlying explanation here is that, under opening electrochemical conditions, the adsorbates engender the metal electrode to undergo a structural transformation, sometimes reversible, of liberating from the substrate and evolving into a more active phase, which includes but is not limited to nanocluster, [127,149,150] transitional surface phase, [151] and undercoordinated structure. [152]In situ EC-STM techniques monitored reversible CO-adsorption-induced local morphological changes of Cu (111) with flat terraces and smooth edges under a reducing potential of −0.85 V SHE to a threadlike Cu nanostructure and small adatom islands when potential increases from −0.85 to −0.45 V SHE . [152]Theoretical microkinetic simulations evidence that the metastable Cu adatom is much more active than Cu (111) surface in CO oxidation and even comparable with Pt (111) (Figure 8e), benefiting from the enhanced electronic back donation from Cu 5d electrons to CO.This structure dynamism is not observed in absence of CO, since the Cu adatom tends to aggregate if not stabilized by CO, which is identified by the very newly released work of Cu (111) electrocatalyzing HER in acidic electrolyte. [127]The strong H adsorption and the resulting high coverage under reductive potential engender the Cu (111) configurational transformation from a pristine metal surface to an organized (4×4) superstructure of Cu adatoms with 12 H adsorbates in the unit cell, combining triangular and square arrangements, which is captured by operando ECSTM at the potential range of −0.25 to −0.27 V versus SHE at pH = 2. Cheng et al. utilize global optimization algorithms and GC-DFT calculations to efficiently identify the active potential window for HER on the Cu (111) surface in an acidic solution lies in the stable potential range of −0.42 V to −0.58 V of the formed low coordination Cu adatoms (A1-15) (Figure 8f). [127]The restructured A1-15 Cu adatoms with H show much better electro-performance with lower reaction barriers toward HER than pure Cu (111) metal surface (Figure 8g).
In acidic media, the adsorbate *H can not only drive the Cu adatoms evolving but also induce dramatic restructuring of the formation of stripes on Cu (100) metal surfaces at ≈−0.4 V, as observed by in situ STM. [153]The enormous chemical space of Cu (100) spanned in structural and configurational dimensions under varying applied potentials is explored efficiently by global optimizations using GCGA, which not only reveals the origin but the structure evolution of the reported stripe pattern. [128]The row-shifting, triggered by H adsorption and negative potential (Figure 8h), is promoted thermodynamically and kinetically as potential decreases and H coverage ( *H ) increases (Figure 8i), originating electronically from the weakening of Cu-Cu bonds between the top-and sub-layer by the H migration by chemical bonding analysis (Figure 8j).The detailed stripe formation pathway with multi-intermediate stages as H coverage varies is uncovered atomically.

The Adaptive Response of Electrolytes and Reaction Species
Apart from the dynamic behavior of the electrode, water molecules are sensitive to the reaction conditions.The dynamic coordination behaviors of the SA site are reported on FeNC and MnNC in the course of ORR. [13,148]Chen et al. carried out AIMD on FeNC with a fully explicit solvation shell to discover solvent water can dynamically coordinate to the Fe center in the backside axial position to adapt the binding strength and coordination configuration as the reaction proceeds, to module the ORR energetics and structural stability of the MNC moiety by varying the Fe-O axial bond lengths (Figure 9a). [13]In particular, the high potential-induced leaching observation is a big concern of the stability of MNCs in long-term operation. [147,154]However, the backside axial water stabilizes the Fe center through the octahedral configuration by alleviating the out-plane distortion of SA sites.Resorting to a simplified gas phase model with or without the backside axial water ligand, a significant geometric difference in the FeN 4 moiety is observed.(Figure 9b) For instance, in the *OH intermediate without backside water, the Fe center is dragged out of the carbon matrix plane by the hydroxyl dramatically to form a square pyramidal configuration, while in the *OH with backside axial water the Fe is relatively kept in-plane in a slightly distorted octahedral configuration.In addition, the distortion angle of Fe, induced by the front-side species, is largely remedied in the presence of the opposite water (Figure 9c).On the other hand, the adaptive behavior of the backside water molecule varies with different adsorbate species as ORR proceeds, engendering the Fe SA featuring distinct oxidation states (Figure 9d).The electronic structure analysis of the Fe center reveals that the predicted Mulliken spin population is S = 2/2 for Fe(IV) in *O…OH − .This emphasized the appearance of the backside water altering the crystal field of Fe to an octahedral pattern since the square planar one would predict S = 0 for Fe(IV).More critically, a novel dissociative ORR mechanism to form "pseudoadsorbed" OH species is discovered unprecedently, as aforementioned (Figure 1c,d).In the PCET steps, the proton species (in the form of hydronium in neutral/acidic media or water in an alkaline medium) can protonate the pseudo-adsorbed hydroxide without traveling to the direct catalyst surface.This, therefore, expands the reactive region beyond the direct catalyst surface, boosting the reaction kinetics via alleviating mass transfer limits (Figure 9e).
The reaction species of OH − in water-electrolyte can work as a modulator as well.A microkinetic model (MKM) is employed to reveal the self-adjusting behavior of FeNC and MnNC. [148]The active center, Fe site, is covered with the intermediate *OH ranging from 0.28 to 1.00 V (Figure 9f).Distinctively, such *OH becomes part of the active center moiety, Fe(OH)NC, optimizing the electronic structure of SA Fe to boost ORR, exhibiting a selfadjusting mechanism and predicting a theoretical half-wave potential of ≈0.88 V.Such intrinsic environmental suitability of SACs demonstrates the necessity of assessing the effect of intrinsic intermediates and the significance of the explicit incorporation of water solvation in the single-atom electrocatalytic process.][157] By combining in situ Raman spectroscopy and AIMD simulations, three evolutional characteristic water configurations in interfacial region, namely, 'parallel', 'one-H-down,' and 'two-H-down', have been determined as the electrode potential decreases at electrified Au single-crystal electrode surfaces, arising to the destruction of the H-bonds network. [34]The statistics of H-bonds rationalizes the observed blueshift of the O-H stretching mode ( O-H ) in Raman spectra (Figure 9g).Subsequently, the open question here is how the physical orientation of water molecules impacts the electrochemical performance of the electrode materials.Apart from the dynamic coordination behavior of interface water on SA sites and Reproduced with permission. [13]Copyright 2022, Springer Nature.f) Coverages of the most abundant ORR intermediates on the SA Fe site of the FeNC center as a function of U. Reproduced with permission. [148]Copyright 2019, American Chemical Society.g) Calculated number of hydrogen-bond donors (N donor ) of interfacial water (red circles) as a function of potential, in comparison to the experimental Raman frequency  O-H (blue circle).Reproduced with permission. [34]Copyright 2019, Springer Nature.
their orientation, the chemisorption of water on metal electrodes is key to differential Helmholtz capacitance. [158]AIMD calculations by Le et al. suggest that at more positive potentials, the interface water will adsorb on the surface of Pt (111), which follows the Frumkin adsorption isotherm.The chemisorbed water can induce a significant interface dipole potential and lead to a potential change and, hence, a negative capacitive response.Combining with the normal dielectric response of the solvent, the experimentally observed bell-shaped differential capacitance of the Helmholtz layer can be reproduced and understood at the molecular level.

Alkali Cation Effects and Interfacial pH
Particularly, the presence of alkali metal cations (M + ) in the interfacial region [159] can significantly alter the electrochemical reactivity and selectivity by directly participating in the reaction [160,161] or by regulating the local electric field, as suggested by Chan group, [39,[162][163][164] or by tuning local pH [40] to modify the kinetics.Conventionally, the concentration of cations is determined by their size, in which smaller sizes lead to higher local concentrations, resulting in a greater accumulated surface charge and the subsequent interfacial field under a certain potential.With  and K).Reproduced with permission. [165]Copyright 2022, Springer Nature.c) Calculated values of cathode pH (left), and cathode CO concentration (right) over Ag electrodes at −1.0 V versus RHE in CO 2 -saturated 0.1 M MHCO 3 (M = Li, Na, K, Rb, Cs) electrolyte.d) Current densities (left) and FEs of CO and H 2 (right) versus applied potential on Ag cathode at −1.0 V versus RHE.Reproduced with permission. [40]Copyright 2016, American Chemical Society.e) Coordinate number change of K + to the *CO 2 and *OCCO, and that in the Bader charge (Bader q) of those two intermediates during DFT-CES iterations at −0.5 V SHE .f) Cation coupling to the intermediates g) A schematic showing the energy level change of CO 2 with and without metal cations, based on the PDOS analyses.Reproduced with permission. [160]Copyright 2022, Springer Nature.h) OCO angle in CO 2 and its Bader charge in the presence of different metal cations.Reproduced with permission. [161]Copyright 2021, Springer Nature.
an exception of the hydroxyl adsorbate on the platinum surface in the process of HER, DFT calculations have evidenced that with the presence of alkali cations, the surface electron redistribution locates between the adsorbates *OH and the nearby water molecules, instead of the Pt surface, and follows the trend of Li + > Na + > K + (Figure 10a). [165]The induced stronger local electric field from its higher charge density in Li + solutions is responsible for the more distinct charge redistribution.Briefly, the C dl is limited by the potential-induced *OH surface coverage and *OH adsorption strength.At a higher potential regime (> 0.60 V RHE ), where there are rich *OH species on Pt, the hydration-sphere size and the related cation-surface distance play a critical role in the capacitance (C dl (K + ) > C dl (Na + ) > C dl (Li + ) while at lower potential region (< 0.60 V RHE ), where there is only a limited *OH coverage, the capacitance (Li + > Na + > K + ) is more dictated by the local cation concentration, echoed with the common sense (Figure 10b).Generally, cations from the electrolyte serve as counterions to balance the charge of the negatively charged electrode. [14,88]However, the direct adsorption of alkali cations on the electrocatalyst of one-dimensional cobalt-dithiolene metal-organic frameworks (MOF) is worth our attention. [166]DFT calculations have demonstrated that the hydration energy of the electronegative Na + is much higher than that of the adsorption energy on the bridge-sites between the two S atoms.The preadsorbed Na + induces the charge redistribution between the S atoms, which has a profound effect on the ΔG *H values.Li with higher electronegativity than that of Na and K has the lowest ΔG *H .
On the other hand, the local interfacial pH is an echo of cation identify (Figure 10c). [40]Experimental investigations show the interfacial pH decreases ≈9 to ≈7 with increasing cation radius from Li + to Cs + over Ag and Cu electrodes, accounting for the buffering capability and the relevant electrochemical CO 2 RR activity of Li + < Na + < K + < Cs + in basic conditions, corroborated Zhang et al. observation of the interfacial pH trend in a bicarbonate electrolyte in CO 2 RR by utilizing a rotating ring disc electrode. [167]However, they claim that the CO 2 RR to CO activity remains unperturbed as local pH varies without further elaboration.As perspectives vary, no unanimous conclusion can be drawn currently.This dilemma in pH perspective by cations calls for a comprehensive investigation of the cation status in the electrochemical process (Figure 10d).
It is dictated that metal cations stabilize certain CO 2 RR intermediates or intimately contact the adsorbates through local electrostatic interactions within the electrical double layer. [39,160,161,163]pecifically, Shin et al. focused on the alkali cations' coordinating ability to all possible intermediates in CO 2 RR on Au (111) and Cu(100) electrodes, empowered by DFT in classical explicit solvent (DFT-CES) methods, and proposed a cation coupled electron transfer (CCET) based mechanism for CO 2 RR, in which K + assists the electron transfer between the electrode and the adsorbed species. [160]The increased coordination number of K + to *CO 2 and *OCCO and the decrease in the intermediates Bader partial charges during iterations suggest a direct coordination stabilization (Figure 10e).Schematic projected density of states (PDOS) shows a downshift of the *CO 2 LUMO after the cation-coupling.Also, the Koper group manifests the obligated mechanistic role of M + in CO 2 RR to CO, making them propose an M + -coupling to *CO 2 − mechanism (Figure 10f). [161]AIMD modeling reveals that the partially coordinated M + with surrounding waters stabilize the *CO 2 − intermediate in regards to OCO angle and Bader charge (Figure 10g).
In addition to cations, the anions in acidic media are in fact not innocent to the electrochemical processes.Luo et al. have provided CVs-based experimental evidence that the stronger interaction of the ClO 4 − with *OH layer leads to slower kinetics of the *OH ↔ *O transition, and a slower ORR rate. [168]Based on the observations, they concluded that the rate of the *O↔*OH process acts as a descriptor of the ORR rate to rationalize the different ORR activity trends for stepped Pt surfaces in acidic and alkaline phases.However, the corresponding interfacial mechanism, in which how anions function the ORR reaction, is still vague and calls for an urgent theoretical understanding to clarify microscopic details.
In this section, we demonstrate with examples that the efficacy of contemporary computational techniques in describing the structural and dynamic features of the solid-liquid interface in electrochemistry.From molecular engineering perspectives, the incorporation of theoretical simulation unveils the identity of the active site of electrocatalysis and constructs a roadmap for structure-activity relations.The underlying electronic configurations of the active sites and the supporting substrate provide a fundamental understanding for the micro-/macro-electrochemical performance.From the interface-modeling viewpoint, AIMD simulations verify the catalytic interfaces should be surveyed as an evolving statistical ensemble of multiphases, with a focus on the electrode surface and SA sites, electrolyte water molecules, and reaction species.Additionally, researchers should pay attention to other intricacies involved in the electrochemical courses, such as the distribution of metal cations.

The Conclusions and Perspectives
Heterogeneous electrocatalysis has emerged as a powerful method in energy conversion processes, with impressive catalytic performance.This review summarizes fundamental concepts in heterogeneous electrocatalysis and detailed thermodynamic and kinetic modeling methods.We also provide a few selected examples from our current research activity, demonstrating how the state-of-the-art AIMD simulations can offer a microscopic understanding of properties at solid-liquid interfaces and effectively guide experiments to aid in describing the structure and dynamics at interfaces, as well as understanding the microscopic origin of the measured reactivity in experiments.Our contributions underline the significance of considering explicit solvation shells in simulating the electrochemical interface, the new proposed descriptor derived from charge-dipole interaction to evaluate HER kinetics on SA sites, and the electrode potential-determined reactivity and selectivity on various electrochemical processes and surfaces.Additionally, we systematically emphasize the unique structural and electronic properties of various surfaces, such as the regulation of the coordination environment of active sites and the corresponding electronic configuration of the active center, aiming to provide essential guidelines for screening and designing electrocatalysts at the atomic scale.We focus on the dynamism of interfaces, including the leaching of SA sites, surface restructuring of metal electrodes under opening conditions, and the response behaviors of the electrolyte molecules to the applied potential, as well as the roles of alkali metals and interfacial pH effects on reaction kinetics.We conclude that theoretical modeling can lead to a new understanding of the electrocatalytic interface and more efficient interfacial design.In addition to the general challenges of making AIMD techniques more accurate and efficient, which are common to other fields, such as solid-state physics or computational biochemistry, we would like to discuss a few specific points that we believe are specific to solid-liquid interfaces.
Despite significant progress in modeling solid-liquid interfaces, several challenges still need to be addressed to develop more realistic and predictive models.One of the primary challenges we foresee is the careful identification of the genuine active center in heterogeneous electrocatalysis.The support phase plays a crucial role in this regard.Theoretical observations suggest that active sites may exist in the adjacent C site of metal Ni in NiN 1 , NiN 2, and NiN 3 electrodes in the Volmer reactions. [88]urthermore, the direct participation of W sites in adsorbing H on Pt ads @WC X alerts us to the impact of the substrate phase in the electrochemical process. [120]It is crucial to develop substrate materials with large specific areas, good wetting abilities, and durable mesoporous structures to enhance the performance and applicability of the active phase.However, equal attention should be paid to the role of supports as to the supported active center at the interfacial region.
Another challenge in modeling solid-liquid interfaces is that the pre-existing mechanisms and descriptors cannot be generalized for all electrode surfaces.We should not apply mechanically to the evolving electrochemical solid-liquid interface.New reaction pathways could be observed when taking comprehensively into consideration of various factors in simulating the interface.Such as the fully explicit model systems help us to identify the "pseudo-adsorbed" hydroxide species in ORR, [13] dual-site participation on NiNC rationalize the observed lower overpotential than that of FeNC and CoNC in OER, [12] the change of local coordination and position of SA sites in MNC helps the formation of a stable dihydride or dihydrogen intermediate (HMH) in HER. [81]In addition, in view of the screening of the electrocatalysts, the descriptors established on bulk metal electrodes could be less effective for SA sites.The traditional volcano plots have been demonstrated to fail to describe the Volmer reaction on MNC due to the strong charge-dipole interaction between *H and with surrounding water molecules. [81,85]The oversimplicity of the capacitor models for 2D materials in the presence of polarizable implicit solvation in conventional CPM, however, arises non-negligible deviations. [52]he stability of the electrocatalysts is a significant concern in the interface region.As discussed in previous sections, the undercoordination and adaptive bonding with the substrate under electrochemical conditions could engender the SA sites or the metallic surface atoms destabilize and eventually aggregate into larger nanostructures.This will not only affect catalytic activity by favorably binding certain reaction intermediates but also compromises the loading efficiency and purity of SAC, hindering widespread commercial applications.This non-negligible issue also poses complexities in the atomic-scale understanding of the interface and limits the rational design of electrocatalysts.Recent grand canonical DFT calculations have demonstrated the existence of an electrochemical potential window within which the aggregated metal target can leach away and the resultant SACs or SCCs are more stable than NPs. [169]Here comes the challenge that seeking the balance conditions between maintaining the stability of the electrocatalysts or avoiding the poisoning of the instantaneous active center by reaction species to gain simplified and unambiguous information, and simultaneously maximizing the catalytic power calls for a more realistic description of the catalytic interface.
Although many experimental and theoretical works have demonstrated the different alkali metal cations at the interface facilitate the electrocatalytic processes by modifying the local electric field, buffering the interfacial pH, or stabilizing the reaction species, [170] a full understanding of this enhancement effect has been a topic of considerable debate.The potassium cations have been manifested to accelerate the CO desorption and CO 2 activation by increasing surface charge density [159] while CO 2 RR is reported to be completely inhibited in the absence of metal cations in solution. [161]Therefore, to more completely understand the descriptor that dictates the electrocatalytic activity, it is essential to uncover how the given cations alter the local chemical environment, i.e., the underlying microscopic structure and evolutional behaviors of the ions at the electrified interface.What's more, to the best of our knowledge, the interfacial structure and electrostatic potential distribution of the condensed cation layer in the EDL are only accessible to state-of-the-art machine learning (ML) due to limitations in spatial and time scales in AIMD, and general design principles are still insufficient.
Regarding computational capability, the combination of DFT and AIMD does help us well interpret the solid-liquid structures and rationally optimize a predictive model.However, the illimitable chemical space of the configurational and dynamical interface, coupled with the infinite time dimension of molecular dynamics is computationally unattainable with AIMD.Recently, the development of ML potentials has significantly revolutionized atomistic modeling in a more extended dimensionality without compromising the ab initio accuracy.Benefitting from systematic ML accelerated MD (MLMD), accurate and efficient calculation of fundamental thermodynamic factors in electrochemistry, including redox potentials, acidity constants, and explicit solvation-free energies have been accomplished by the designed automated workflow with several orders of magnitude speedup. [171]In this regard, as ML continues to develop, it is bound to play an increasingly indispensable role in electrochemistry and other fields.[174][175] The collaboration of theoretical and experimental techniques is an invincible approach for modeling electrocatalysis.However, it is time to retrospect the connections between the experimental measurements and computed parameters, such as the efficiency of the microkinetic model (MKM) to relate inherent energetics of the active sites and experimentally measured electrocatalytic rates, [88] to gain reliable interface insights. [176,177]What's more, in electrochemistry, the transient capture of the structural evolution during operation highlights the future theoretical direction, such as ML, to interpret complex in situ spectroscopic signals [147] while computational high-throughput sampling of electrocatalytic materials requires a sufficiently large experimental dataset to produce validated predictions. [178]s we contribute practical models and methods to probe a more realistic solid-liquid interface, concurrent challenges are confronting us.Computational electrochemistry holds great potential for advancing our atomistic understanding of electrochemical systems and developing new strategies and materials to address pressing energy and environmental challenges.In this review, we have presented current theoretical endeavors, opening challenges, and promising future directions in exploring the solid-liquid interface in electrochemistry.We anticipate this thematic article would inspire new research enthusiasm.

Figure 1 .
Figure 1.a) Selected snapshots from AIMD trajectories during the process of configurational transformation from end-on *O 2 to side-on pattern with liquid water phase on FeNC.b) O-O bond length evolution during the AIMD simulation with explicit solvation.Reproduced with permission.[57]Copyright 2020, American Chemical Society.c) Selected snapshots from AIMD trajectories show an unconventional dissociative ORR pathway that produces pseudo-adsorbed hydroxide species.d) Probability density distribution of coordination number of O by H (CN O ) in Z-direction, from the three equilibrated AIMD trajectories of *O…OH − , *OH…OH − , and *…OH − , showing the spatial distribution of pseudo-adsorbed OH.Reproduced with permission.[13]Copyright 2022, Springer Nature.e) Statistic numbers of H-bonds between the solvation water and CO 2 reactant along the adsorption reaction coordinate at different potentials on FeNC and snapshots of the solvation environment around CO 2 before and at TS. f) HOMO and g) spin density distribution of the CO 2− anion in a linear or bent configuration.Reproduced with permission.[14]Copyright 2022, American Chemical Society.h) Evolution of the O-H distances during proton transfer and the snapshots of the AIMD trajectories for the *COO s H to *COO a H species on Cu metal surfaces.i) Reaction free energy barrier comparison of CO formation form *COO s H and *COO a H. Reproduced with permission.[58]Copyright 2022, American Chemical Society.j) Charge density difference and k) pDOS and pCOHP comparisons of N 2 adsorption with (lower panel) and without (upper panel) water shell.Reproduced with permission.[59]Copyright 2022, American Chemical Society.

Figure 2 .
Figure 2. a) Adsorption-free energies of *OOH and *O as a function of *OH on all studied MNCs.Theoretical and corresponding experimental onset potentials for b) ORRand c) HER versus the descriptor .Reproduced with permission.[67]Copyright 2018, Springer Nature.d) The comparison results of adsorption-free energies of H between PBE+U and the descriptor  for MNCs.Reproduced under the terms of a CC-BY license.[70]Copyright 2022, The authors, published by ACS.e) Relationship between the reaction Gibbs free energy of the CO, CHO, OH, and H 2 O on TM 1 /Cu(111) as a function of *OH binding energy.Reproduced with permission.[66]Copyright 2022, Tsinghua University Press.f) The linear relationship of the first and the last hydrogenation step of NRR on a series of dual-atom catalysts by DFT calculations.Reproduced with permission.[75]Copyright 2021, American Chemical Society.

Figure 3 .
Figure 3. a) RDF analysis of directed H 2 O molecules around the *H in the sphere of *H•••H distance on different MNCs.b) The illustration of the special *H…HOH moiety.c) The averaged Bader charges of *H on MNCs during AIMD simulations.Analysis of the angle between adsorbed hydrogen and O-H bond of water (*H•••O-H) along the reaction coordinates of H adsorption on d) CrNC, and e) RhNC.Net charge evolution of f) graphene matrix and g) TM sites on different MNCs.h) Theoretical evaluation of the acidic HER activity on different MNCs.Reproduced with permission.[85]Copyright 2023, American Chemical Society.

Figure 4 .
Figure 4. a) Pourbaix diagram showing the thermodynamic equilibria of the CPET, ET, and PT steps scale differently with pH.Reproduced with permission.[91]Copyright 2016, the Royal Society of Chemistry.b) Rate map at −0.8 V SHE and pH = 2 for CO 2 RR to CO across electrodes.Reproduced with permission.[93]Copyright 2021, Springer Nature.c) The competing relation of the formation of H 2 O 2 and H 2 O in the ORR process on CoNC in perspectives of thermodynamic (upper panel) and kinetic (lower panel) modeling.Reproduced with permission.[96]Copyright 2021, American Chemical Society.d) Potential-dependent free energy profiles of CO 2 adsorption at the FeNC-water interface.e) Location of transition states (TS) during adsorption at different potentials.The fitting linear relationship between f) ΔG and ΔG ‡ and g) ΔG/ΔG ‡ versus potential.Reproduced with permission.[14]Copyright 2022, American Chemical Society.h) Potential-dependent free energy profiles of CO 2 adsorption (upper) and *COOH formation under basic (lower left panel) and acidic (lower right panel) environment on NiNC.i) LPD-K predicted partial current densities for CO evolution in the basic (left) and acidic (right) phase.Reproduced with permission.[88]Copyright 2022, American Chemical Society.j) Free energy profiles for *N 2 hydrogenation to *NNH on FeNC.Reproduced with permission.[59]Copyright 2022, American Chemical Society.Bader charge evolution in CO 2 chemisorbing process on k) NiNC.Reproduced with permission.[88]Copyright 2022, American Chemical Society.and l) FeNC.Reproduced with permission.[14]Copyright 2022, American Chemical Society.
reaction barriers of the key elementary steps of *OOH → *O + OH − and *OOH + H 2 O → * + H 2 O 2 + OH − , which demonstrated that the *−O bond breaking has a lower barrier than the *O−OH breaking.However, the thermodynamic information based on the CHE model gives a completely different conclusion that the strong thermodynamic preference for the O−OH bond breaking and the H 2 O formation (Figure

Figure 6 .
Figure 6.a) The theoretical structures of Pt doped @WC x and Pt ads @WC x .b) Calculated Gibbs free energy diagram for water dissociation on different electrocatalysts.c) Projected densities of states (PDOS) of Pt 5d orbitals of Pt(111), Pt doped @WC x , and Pt ads @WC x .d) Gibbs free energy of *H adsorption for different catalysts.e) The Bader charge analysis of the atomic Pt site.f) PDOS of Pt 5d orbitals and W 5d orbitals in the Pt doped @WC x and Pt ads @WC x .Reproduced with permission.[120]Copyright 2022, Wiley-VCH GmbH.g) Free energy diagram for ORR reaction pathways on Cu 1+ -SA and Cu 2+ -SA active sites in alkaline media.(pH = 13) h) PDOS overlap for Cu(I) and Cu(II) with *O atoms.Reproduced with permission.[121]Copyright 2020, American Chemical Society.i) Linear correlation of the CO 2 activation barrier, ΔG ‡ , versus U RHE on NiNx electrodes.j) Bader charge of Ni center for NiNx under PZC.k) PDOS of 3d x2-y2 orbitals for Ni center in all sites and C 2* orbitals of CO 2 .Reproduced with permission.[88]Copyright 2022, American Chemical Society.

Figure 7 .
Figure 7. a) Schematic illustration of orbital interactions between adsorbed CO (5 and 2*) and 3d orbital (d z2 , d xz /d yz ) of Fe site in Fe-SAC (upper)and NiFe-DASC (lower).b) Calculated free energy diagrams for CO 2 RR and OER processes.Reproduced with permission.[29]Copyright 2021, Springer Nature.c) Proposed reaction scheme of the single-site (upper) and dual-site (lower) mechanisms towards OER.d) Calculated free energy diagrams at 1.23 V for OER over FeNC, CoNC, and NiNC.Reproduced with permission.[12]Copyright 2018, Springer Nature.e) Calculated free energy diagram of NiPc-CN-H 2 , NiPc-H 2, and NiPc-OMe-H 2 electrocatalysing CO 2 RR at −0.11 V SHE .Reproduced with permission.[124]Copyright 2017, Springer Nature.f) Electron density difference plot showing charge transfer from graphene substrate to NiPc molecules.Reproduced with permission.[126]Copyright 2021, American Chemical Society.

Figure 8 .
Figure 8. a) Illustration of the dynamic reversible transformation between the CuNC and Cu 4 clusters under working conditions.b) Calculated free energy (ΔG H ) of the 1 st , 2 nd , 3 rd , and 4 th H adsorption under different potentials.Cu-N bond lengths evolution for c) the complete SA leachingand d) the incomplete one.Reproduced with permission.[117]Copyright 2022, American Chemical Society.e) Calculated current densities for the Pt(111), Cu(111), and the undercoordinated Cu adatom active site motif generated at the restructuring Cu electrode.Reproduced with permission.[152]Copyright 2020, Springer Nature.f) Calculated surface Pourbaix diagram showing potential and pH dependence of the Cu (111) surface state.g) Reactive evolution of the A1-15 structure.Reproduced with permission.[127]Copyright 2023, Wiley-VCH GmbH.h) Calculated surface Pourbaix diagram showing potential and pH dependence of the Cu (100) surface state.i) Total grand canonical free energy as a function of H coverage at different electrode potentials.j) Color-filled contour map of Cu(100) with 1 H showing a weakening of the Cu-Cu bond between the top layer and the sublayer.Reproduced with permission.[128]Copyright 2022, American Chemical Society.

Figure 9 .
Figure 9. a) The radial distribution function g(r) of the backside axial water on FeNC in different reaction intermediates with representative snapshots labeled on the top.b) Overlapped geometries of gas phase reaction intermediate with (colored) or without (transparent blue) backside axial water ligand.c) The evolution of FeNC distortion angle along the reaction coordinate for gas phase reaction intermediates with (blue lines) or without (orange lines) backside water ligand.The shaded area represents the difference of FeNC distortion angle in the two cases.d) The scatter plot of Mulliken spin population versus the Bader charges on Fe along the reaction coordinate.e) The z coordinate of the O in dissociated OH − relative to that of the Fe single site (denoted as R) in *O…OH − , *OH…OH − , and *…OH − intermediates.Reproduced with permission.[13]Copyright 2022, Springer Nature.f) Coverages of the most abundant ORR intermediates on the SA Fe site of the FeNC center as a function of U. Reproduced with permission.[148]Copyright 2019, American Chemical Society.g) Calculated number of hydrogen-bond donors (N donor ) of interfacial water (red circles) as a function of potential, in comparison to the experimental Raman frequency  O-H (blue circle).Reproduced with permission.[34]Copyright 2019, Springer Nature.

Figure 10 .
Figure 10.a) Electron density difference maps of the Pt (111)-OH ad -water interface after introducing Li + (left), Na + (middle), and K + (right).b) Experimental double-layer capacitance (C dl ) on the Pt disc electrode at different applied potentials in 0.1 M MOH solutions (M = Li, Na,and K).Reproduced with permission.[165]Copyright 2022, Springer Nature.c) Calculated values of cathode pH (left), and cathode CO concentration (right) over Ag electrodes at −1.0 V versus RHE in CO 2 -saturated 0.1 M MHCO 3 (M = Li, Na, K, Rb, Cs) electrolyte.d) Current densities (left) and FEs of CO and H 2 (right) versus applied potential on Ag cathode at −1.0 V versus RHE.Reproduced with permission.[40]Copyright 2016, American Chemical Society.e) Coordinate number change of K + to the *CO 2 and *OCCO, and that in the Bader charge (Bader q) of those two intermediates during DFT-CES iterations at −0.5 V SHE .f) Cation coupling to the intermediates g) A schematic showing the energy level change of CO 2 with and without metal cations, based on the PDOS analyses.Reproduced with permission.[160]Copyright 2022, Springer Nature.h) OCO angle in CO 2 and its Bader charge in the presence of different metal cations.Reproduced with permission.[161]Copyright 2021, Springer Nature.