All-atom XL-BOMD simulations show oxygen reduction on nitrogen-doped graphene switching between inner- and outer-sphere mechanisms with applied bias

The work provides the first demonstration of XL-BOMD applied to heterogeneous electrocatalysis, simulating electrode and electrolyte together with potential biasing in a unified framework, running stably up to 20 ps with a 0.1 fs time step and about 0.15 s wall-clock per MD step for a 622-atom, 1840-electron NG/water/O2 system. Prior static DFT with the computational hydrogen electrode assumes proton-coupled electron transfer and cannot describe non-adsorptive charge transfer, while conventional quantum molecular dynamics is limited by the cubic scaling of Kohn-Sham DFT and self-consistent field convergence; XL-BOMD treats electronic degrees of freedom as dynamical variables and uses a shadow potential to avoid iterative SCF at every step. The paper reports no energy drift over 10 ps NVE simulations, timing scaling of GPMDK with system size, and second-order power-law SCF convergence with preconditioning; DFTB adsorption energies, water conformation, and the O2 bond length (1.24 Å) are compared against DFT and experiment.

The bias applied to the NG sheet determines whether the ORR proceeds by an inner- or outer-sphere pathway: at higher biases (0 to 3.0 eV) O2 is reduced without chemisorption via solvent-mediated electron transfer, with four electrons transferred from the NG sheet in two steps and four OH- formed; at lower biases (-1.25 to -0.25 eV) O2 adsorbs on the NG surface before direct electron transfer. The paper introduces a method to apply an electrochemical bias by shifting the onsite orbital energies of the NG atoms, distinct from constant-potential methods, keeping the ensemble canonical (NVT) so the reaction driving force can be scanned systematically. Conclusions rest on 58 production simulations up to 20 ps, covering 10 inner-sphere and 20 outer-sphere applied biases with three to five replicates each; species and bond cleavage are identified from Mulliken charges and O-O distances, and three control systems (NG in vacuum, O2 in water, pure graphene in water) showed no ORR.

The simulated kinetics are not fully captured by Butler-Volmer or Marcus-Hush-Chidsey theory: both inner- and outer-sphere data can be fit with two Tafel slopes, but no corresponding change in rate-determining step or mechanism is seen, and the Tafel slope bends continuously; MHC theory fits the full dataset convincingly even though the simulations clearly show two distinct mechanisms. This cautions that inferring mechanism from macroscopic kinetic fits alone can be misleading and that an atomistic, step-by-step picture is needed to complement such theories. Evidence comes from reaction-time statistics over 58 simulations (reaction time defined as the time to form all OH- ions) and fitted slopes of 49.2 and 189 mV dec-1 for inner-sphere and 49.3 and 316 mV dec-1 for outer-sphere; the paper also notes the ranges do not span the one decade of current typically preferred for Tafel analysis.

The electronic-structure origin of NG's higher activity is made concrete: relative to pure graphene, nitrogen doping both shifts the electron chemical potential to increase the reduction potential and increases the density of states around the Fermi level, providing a source of electrons; after ORR the NG sheet is oxidized, creating four electron holes mostly localized on carbons directly bonded to nitrogen and delocalizable across the sheet via resonance structures. This links the macroscopic catalytic advantage of NG over pure graphene to computable density-of-states features and hole distributions. Based on density-of-states comparisons among O2, carbons near nitrogen in NG, and pure graphene, and on average Mulliken charges of neighboring carbons over 0.5 ps after oxidation (first-neighbor carbons about 0.202, second-neighbor about -0.019).

This work is aimed at researchers in theoretical electrocatalysis, computational electrochemistry, and quantum molecular dynamics methods, and applies to settings where potential-dependent reaction mechanisms at solid-liquid interfaces are studied with explicit solvent. It enables researchers to treat electrode and electrolyte together in one simulation, apply a controlled bias, and observe electron-transfer events step by step, and the authors propose extending it to steady-state current conditions, explicit pH effects, and defective or layered catalyst architectures. The results concern isolated single-molecule ORR events on an ideally polarizable electrode, deliberately avoiding coverage effects, adsorbate-adsorbate interactions, and competitive multi-site binding.

Several open questions remain for a careful reader: the bias-to-overpotential conversion is stated to need further testing and validation; the simulation conditions are described as most likely corresponding to high overpotentials and use acceleration strategies such as placing O2 near the surface initially and steered molecular dynamics; the final step of the inner-sphere pathway (adsorbed oxygen reacting with water to form two additional OH-) was not observed, which the authors hypothesize relates to partial reduction or simulation length; transient species assignments rely on Mulliken charges compared with reference tables, and the authors note Mulliken charges do not formally indicate oxidation state and that adsorbed O2- and O*- species have similar charges in isolation and cannot be distinguished; and transient species formation in the outer-sphere mechanism may relate to simulation conditions with no steady-state current, no electron reservoir, and no constant bias. In addition, this is a fast parse of the text, so specific numerical details in figures and supporting information should be checked against the original.