Exploring non-adiabaticity to CO reduction reaction through ab initio molecular dynamics simulation

First-principles methods, such as density functional theory (DFT), have been widely used in a variety of electrocatalytic reactions, such as water splitting including oxygen evolution^1–5 and hydrogen evolution reactions,^6–9 CO₂ reduction,^10–13 and solar fuel cells.^14–16 By utilizing the state-of-the-art computational techniques, such as computational hydrogen electrode (CHE) model,¹⁷ nudged elastic bands (NEBs),¹⁸ and Born–Oppenheimer molecular dynamics (BO-MD),¹² DFT calculations enable a detailed free-energy determination of possible reaction paths, including intermediate states, transition states, and reaction barriers. Various characteristics, for example, pH,¹⁹ electrode potential,^10,17 cation induced electric field,^20,21 and electrolyte,^12,14,22,23 have been explored systematically to illustrate the reaction mechanisms and to design high-performance catalysis. Most of these methods, particularly to determine the reaction barriers, are based on the ground-state DFT. By assuming a much slower reaction process compared to the timescale of electronic thermalization, electronic excitations during the related reactions have been ignored; hence, the adiabatic Born–Oppenheimer assumption is used in many studies on electrocatalysis. The non-adiabatic (NA) effect results from the electronic evolution with a finite timescale. However, the NA effect on electrochemical reactions is rarely explored, and we have very limited knowledge on whether non-adiabaticity will affect electrochemical reactions and to what extent it can contribute.

The electronic NA effect is defined as the change of the electronic ground state to excited eigen-states due to the time evolution of the wave functions. This effect manifests itself in a chemical reaction when the timescale of its dynamics is faster to that of the carrier’s (electrons and holes) thermalization. As a result, the excited electron (and hole) is not always at its equilibrium ground state. Meanwhile, for electron transfer related reactions, the NA effect can result in charge transfer bottleneck, where carriers need finite time to move from the carrier donor to the acceptor. In comparison, when the reaction is carried out very slowly, the fast relaxation of electronic excitation or fast charge transfer makes the NA effect insignificant to the reaction. Owing to the ultrafast nature, the NA effect plays a significant role in chemical reactions such as photochemistry,^24,25 collision,²⁶ atom-stopping,^26–28 and electron transfer processes.^25,29 In particular, a great deal of research has been focused on gas-phase catalytic reactions on surfaces to understand the contribution of the NA effect.^30–35 For example, the NA simulation with the fewest-switches surface-hopping algorithm has shown a strong NA effect for spin flipping and transition during the O₂ dissociative adsorption on Al and Pd surfaces, with the estimated rate consistent to the experiments.^36,37 In another study, the reverse process, associative desorption of N₂ on Ru(001), further showed the NA effect from an ab initio simulation.³⁸ In that work, a consistent agreement between the simulation and the experiments could be obtained only after including the NA effect in the calculation. A comprehensive theoretical study has been performed to explore the NA effect of H₂/Cu (110) and N₂/W (110). However, their simulation has shown a marginal effect of non-adiabaticity on the adsorption process of diatomic molecules.³⁷ Based on these examples, the role of the NA effect seems to depend on specific reaction types. However, for electrochemical reactions under aqueous conditions such as heterogeneous catalysis, the NA effect has been rarely studied. Electrochemical reactions necessitate the transfer of charge from one place to another; thus, it is more likely a NA process. Besides the question of an excited state induced by the reaction dynamics, another possibility is the charge transfer bottleneck, which also makes the process non-adiabatic. A recent work focusing on the initial CO₂ adsorption to various metal surfaces has shown very fast electronic hybridization compared to the adsorption,²¹ exhibiting the adiabatic nature of the chemical adsorption process. Different from the initial adsorption process for CO₂ reduction, its subsequent steps involve the proton-assisted electron transfer (PAET). Such a PAET process also exists in water splitting and hydrogen evolution reactions. Some of the fast PAETs can form the transition state within 1 ps or faster,^39–41 comparable to the timescale of the electronic thermalization, indicating the potential role of the NA effect. Nonetheless, the NA effect on electrochemical reactions with the addition of protons is rarely investigated in detail. It leaves many questions unanswered: for example, will the NA effect play a role in electrochemical reactions such as CO₂ reduction involving the fast proton motion? If this is the case, to what extent the change of the reaction barrier can be due to non-adiabaticity? Will the NA effect occur mostly in the carrier excitation or in the charge transfer? How do the electrolyte and the applied electrode potential influence the NA effect?

In this work, we seek to understand the role of the NA effect on heterogeneous catalytic reactions in aqueous conditions. By using one step of the CO₂ reduction ⋆CO + H₃O⁺ + e⁻ → ⋆COH + H₂O (⋆: copper surface) on the copper (111) surface as an example, we perform adiabatic (ground-state BO-MD) and non-adiabatic [Ehrenfest real-time time-dependent DFT (rt-TDDFT)] simulations to study the reaction. In BO-MD, the adiabatic eigen-states are solved by diagonalizing the Hamiltonian at every MD step. The occupation of electrons on each state is based on their eigen-energies. Thus, the electronic structure is always at the ground state for every step. At the contrary, rt-TDDFT evolves as the time-dependent wave function following Schrödinger’s equation, allowing the electronic structure to be excited. Meanwhile, the excitation of the electronic structure might drive the movement of ions differently compared to the ground-state electronic structure. Hence, rt-TDDFT has been widely used to simulate various NA processes such as optical excitations,⁴² proton-assisted chemical reactions,⁴³ and ion sputtering.^27,44 Different from other TDDFT algorithms, where a very small time step (Δt ∼ 0.001 fs) has to be used to evolve the charge density, the implementation we have adopted here uses the adiabatic states [ϕ_j(t)] as the basis to expand the electronic wave function [ψ(t)]. These adiabatic states are solved from the Hamiltonian at each ion’s step t_n with the time step Δt = t_n+1 − t_n ≤ 0.1 fs. The time-dependent wave function is expanded as follows:

where the adiabatic state ϕ_j(t) is solved by diagonalizing the Hamiltonian H at step t: H(t)ϕ_j(t) = ϵ_j(t)ϕ_j(t). The coefficient C_ij(t) for the wave function is evolved from t_n to t_n+1 non-adiabatically, following the Schrödinger’s equation i∂ψ(t)/∂t = H(t)ψ(t) using a much smaller time step. The Hamiltonian in the Schrödinger’s equation is based on the linear interpolation from H(t_n) to H(t_n+1). However, since the adiabatic states (with about 100 states) are used as a basis to construct the wave function and Hamiltonian (instead of plane-waves), the evolution of the wave function from t_n to t_n+1 involves only a small size matrix, and its cost becomes negligible. This method allows us to evolve wave functions and ions’ dynamics of a complex system with hundreds of atoms such as the surface chemical reaction presented here. In this work, the reaction is simulated with a CO molecule adsorbed on the copper (111) surface, and it is attacked by H₃O⁺ to form COH on copper. We find that the BO-MD and rt-TDDFT with the same initial setups could reveal opposite results: near the reaction barrier, BO-MD allows the reaction to happen, while rt-TDDFT fails to proceed the reaction to form ⋆COH but returns back to ⋆CO. Such a difference clearly demonstrates the role of non-adiabaticity on electrochemical reactions in aqueous conditions. The reaction barrier change caused by non-adiabaticity is estimated, which is up to 6% compared to the ground-state method barriers. We also explored the reaction with the explicit solvent model and with different electrode potentials, and we found similar NA effects.

The PWmat^45,46 package based on the plane-wave pseudopotential DFT is used to perform the total energy calculation, structural optimization, BO-MD, and rt-TDDFT. SG15 pseudopotentials⁴⁷ with 50 Ryd plane-wave energy cutoff are used to ensure the convergence of charge density. For structural relaxation, the residual forces are set below 0.02 eV/Å. In our calculation, we choose a four-layer Cu(111) slab surface with a 3 × 3 in x-y direction superlattice. The most widely used DFT exchange-correlation functionals such as PBE and LDA predict CO adsorption on the hallow site of three Cu atoms, contradicting the experimentally observed top-site (on the top of one Cu atom).^48,49 Instead, we choose revised PBE (revPBE)⁵⁰ to reproduce the correct adsorption site. In the dynamical simulations, the reaction barrier is sensitive to the length of the time step, and we find that 0.1 fs is enough to obtain an accurate reaction barrier (see the supplementary material). The rt-TDDFT calculation is also converged with this time step. Both spin-polarized and spin-unpolarized calculations are performed, but these two types of calculations yield quite similar results.

The solvent has been shown to play a significant role in CO₂ reduction reactions on metal surfaces. Two types of solvent models are commonly used to represent the solvent effect in DFT calculations: the implicit solvent model with a continuum dielectric response and the explicit solvent model with water molecules in simulation. In our calculation, we have tried both methods to examine the effect of non-adiabaticity. For the implicit solvent model, when the system contains charged ions such as hydronium with strong solvation energy, the solvent model has to be tuned carefully to yield a correct energy for the ions, so that the energetics of the transition from free hydronium to ⋆COH will be correct. Here, the continuum polarizable solvent model is used with specific ion–solvent interaction parameters.^22,23 The solvent parameters of H and O (belonging to hydronium) are tuned to reproduce the solvation energy of the charged ion. However, computing solvation energy of a charged ion in explicit water is a non-trivial task due to water fluctuations. Instead, we borrow the idea of the computational hydrogen electrode to compute the free energy of the ion under aqueous conditions.⁵¹ Using hydronium as an example, the equilibrium reaction $H3O+(aq)+e−→12H2(g)+H2O(aq)$ happens at potential U = 0 V. Thus, the enthalpy of H₃O⁺(aq) can be expressed as H(H₃O⁺(aq)) = 1/2H(H₂(g)) + H(H₂O(aq)) + 4.44 eV. Here, H stands for enthalpy, and H(H₂O(aq)) − H(H₂O(g)) = −0.274 eV. The water solvation energy obtained from the experiment is 0.274 eV, which has been widely used to fit the solvent model of water.^22,52,53 The hydrogen electrode potential in reference to vacuum is 4.44 eV. Note that the explicit solvent model is used only to describe the enthalpy, instead of the free energy of H₃O⁺(aq), in agreement with the earlier work of implicit solvent model development.^51,52,54 We tune the solvent parameters of H and O, so that the DFT calculated energy of the hydronium ion with the implicit solvent model matches H(H₃O⁺(aq)) obtained with the aforementioned formula.

Figure 1 shows the optimized initial and final structures. The initial structure is built with hydronium close to ⋆CO with a hydrogen bond, which is a local minimum structure. Such a hydrogen bond is optimized, yielding a bond length of about 1.6 Å. To simulate the reaction with PAET using MD, an initial velocity is added to the hydrogen atom of the hydronium close to ⋆CO, with the direction of velocity pointing to the oxygen of ⋆CO. By tuning the magnitude of the initial velocity, we can monitor when the proton can transfer from hydronium to ⋆CO instead of returning. Such an initial kinetic energy of the proton can be treated as the reaction barrier. More precisely, to find the initial atomic configuration and velocity for this reaction to happen at the exact required kinetic energy, we have adopted a “reversed process” procedure. In this procedure, the nudge elastic band (NEB) calculation is performed first to reveal the reaction path and transition state. Then, by starting from the transition-state structure with a very small initial velocities perturbation toward the initial reaction direction, a short BO-MD is performed. This will yield an initial atomic structure. Starting from this atomic position, with reversed velocities, the BO-MD will drive the system to the transition state due to time inversion symmetry. Thus, a slight increase in the initial velocity can lead to a transition to the final state. On the contrary, a slight reduction (e.g., 0.1%) in velocity will prevent the reaction from happening. Using this way, we can quickly identify the adiabatic reaction barrier using BO-MD. For the reaction ⋆CO + H₃O⁺ + e⁻ → ⋆COH + H₂O, the energy difference ΔE = E_final − E_initial is calculated to be around 0.5 eV. In the experiment, an overpotential is added to overcome ΔE or to make it negative to make the reaction to proceed spontaneous. To mimic the applied overpotential to the electrode, we add two electrons to the system and relax the structures so that the energy difference between the initial and final structures is close to zero. Figure 2 shows the reaction paths computed with BO-MD and rt-TDDFT. In this figure, both calculations of BO-MD and rt-TDDFT start from the same initial structures and velocities as well as the initial electronic structure. The initial kinetic energy of the proton equals the BO-MD reaction barrier to just let the reaction happen. For the BO-MD simulation, the proton of hydronium move from H₃O⁺ to ⋆CO quickly at the beginning. Then, it starts to slow down from 10 fs to 25 fs. Eventually, it bonds to ⋆CO after around 30 fs indicated by the exchange of the distances toward CO and H₂O. We extend the simulation up to 50 fs to make sure the proton will not return back to the water molecule. However, rt-TDDFT reveals a completely different reaction path. The proton follows almost the same reaction path of BO-MD at the beginning. But, it deviates from the BO-MD’s path after around 5 fs, proceeding to the opposite results in the end. During the simulation, the proton does not move across the reaction barrier, but it returns back to water molecules re-forming the initial structure. We also perform rt-TDDFT up to 50 fs to confirm that the reaction does not happen during this time.

To unveil the underlying reason for the dramatic difference, Fig. 3(a) shows the eigen-energies for the states near the Fermi level in the BO-MD simulation. The colors indicate the occupation of the states during the reaction. Near the Fermi level, there are two eigen-states with wave functions mostly on adsorbed CO on copper as shown in Fig. 3(b). At time t = 0 fs, these two states are almost degenerate, except that they are split owing to the weak hybridization with H₃O⁺. During the reaction [Fig. 3(a)], most of the states have relatively small changes, except the state hybridized with hydronium near the Fermi level. When the proton is moving close to ⋆CO, the energy of state I becomes lower, indicating the hybridization developed between the proton and ⋆CO. More importantly, we also track the change of the occupation of this state as shown in Fig. 3(d). Initially, at t = 0 and room temperature, state I is 72% occupied, while state II is 55% occupied. As the reaction goes, the occupation of state I rises until it is fully occupied. On the other hand, the occupation of state II slightly reduces. The total occupation of 2.55 electrons on these two states increases to about 3.0 (non-spin case). Thus, there is an increase of around 0.45 electrons on these two levels. Such an electron increase indicates that the previously empty proton is now occupied by electrons after the reaction. This increase in electrons on H manifests the bond formation between H and CO. The major part of the 0.45 electron transfer is provided from the Cu slab. Such an electron transfer from Cu can be verified by a direct charge measurement before and after the reaction. With a horizontal plane (x-y plane) with its z-value in the middle between the top-layer Cu and C atoms, the total electrons above this plane is found to increase by 0.35 after the reaction. This is also consistent with the results reported in Ref. 21. It is interesting that this charge is not 1. Under the computational hydrogen electrode (CHE) approximation, this charge transfer should be 1. It is well known from a fixed total charge canonical calculation, the charge flow from the bulk to the surface layer is less than 1. The reason for less than 1 charge flow is due to screening of the metal surface.

The abovementioned picture is dramatically different in the rt-TDDFT simulation. As shown in Fig. 3(c), at the beginning of the reaction, states I and II change in a similar way as in the BO-MD. But after 15 fs, they become different. More dramatically, the occupations projected on the adiabatic states I and II almost do not change during the simulation time. The occupation on the adiabatic states [f_j(t)] are calculated as $fj(t)=∑iϕj(t)ψi(t)2O(i)=∑iCi,j(t)2O(i)$⁠, where O(i) is the occupation of the time evolution wave function ψ_i(t), which does not change under rt-TDDFT. The charge on H is controlled by both the hybridization strength of the adiabatic CO–H state and the occupation for this state. If the same initial structure and velocities as for BO-MD are used, the relative constant f_j(t) for states I and II by rt-TDDFT leads to less charge transfer from Cu to H, which suppresses the proton’s motion toward CO and reduces the bond strength of the potential CO–H bond. As a result, there is no formation of a CO–H bond (due to the lack of electrons), and the system bounces back to H₃O⁺ as shown in Fig. 2. The lack of charge transfer is also verified by the direct charge measurement above the horizontal plane as discussed earlier. The change of charge from Cu is less than 0.13 electrons, much smaller than that of BO-MD [Fig. 4(b)]. This example clearly shows how non-adiabaticity plays a role in electrochemical reactions. In these calculations, a smearing of 0.1 eV with the Gaussian smearing method is used. By using a room temperature smearing (0.025 eV) with Fermi–Dirac distribution, it yields a similar occupation (particularly for the difference between states I and II), still allowing for NA contribution.

Having seen the aforementioned difference between BO-MD and rt-TDDFT, we intend to see whether the reaction can happen in rt-TDDFT when a larger initial H velocity is provided. Surprisingly, only a small increase in the H initial kinetic energy (12 meV, around 6.1% increase in the corresponding potential barrier) is needed to make the reaction happen. Nevertheless, the underlying microscopic mechanisms might be different. First, the occupations of I and II states are still constant, almost the same as in the case of non-reaction rt-TDDFT simulation with the same initial velocity to BO-MD. On the other hand, as shown in Fig. 4(b), the charge transfer measured from the spatial location using the z-value horizontal plane is almost the same as that in the BO-MD. To solve the puzzle, we have looked at the H₂O–H and ⋆CO–H distance [shown in Fig. 4(a)]. Compared to BO-MD, the ⋆CO–H distance in the reaction of rt-TDDFT case is slightly closer. This small distance, perhaps together with other factors, can cause the change in the adiabatic state. This is shown in Fig. 4(c), where the difference in the charge density of state I [$ψTDDFT2(r)−ψBO−MD2(r)$] is shown as the isosurface plot. We can see that there are more charges in $ψTDDFT2(r)$⁠, above the value-z horizontal plane than $ψBO−MD2(r)$⁠. Thus, although the occupation on state I does not change in rt-TDDFT, the $ψTDDFT2(r)$ itself has changed, inducing the same charge transfer, which is necessary to form the ⋆CO–H bond. We, thus, conclude that it is not necessary to change the state occupation in order to have the reaction happening. We also note that the different state occupations can have similar spatial charge distributions, which also happens in other cases, such as during the d-state valence charge change in transition metal oxides.⁵⁵ 

Finally, we examine the situation with spin-polarization. After turning on the spin, the reaction path shows a negligible difference compared to Fig. 2. We do note that, if the k-point is not sufficient, in some cases, for the BO-MD simulation, after the proton exchange, the system can become spin-polarized. We expect this could be a real case if CO is placed in a small Cu cluster instead of bulk Cu (see Fig. 4 of the supplementary material ). This spin-polarization, however, will never be developed in rt-TDDFT, since such a spin flip is impossible without spin–orbit coupling. Even with spin–orbit coupling, the time of the reaction discussed here will not be enough to make such a spin flip. This indicates that, perhaps, in some single site catalytic reactions (such as a transition metal in the CN monolayer), the NA effect could be bigger.

As discussed earlier, the implicit solvent model reproduces the energetics of the solvation effect to ions. However, it does has its disadvantages:⁵⁶ primarily as an averaged continuum media, it lacks the atomistic bonding information. More importantly, for dynamical simulations, the implicit solvent model has instant dielectric screening response. But, in reality, the surrounding water may not have enough time to rotate itself and rearrange the structure, following the fast proton transfer movement. Meanwhile, the surrounding water molecules could form hydrogen bonds with hydronium or even with ⋆CO to change the energy levels. These hydrogen bonds depend strongly on the detailed structure of the surrounding solvent molecules, instead of an averaged screening potential in the case of implicit solvent molecules. To overcome this challenge, we utilize a hybrid solvent model by sampling an explicit water molecules layer around the reaction site. The implicit solvent model is still used outside the explicit solvent model layer. Obtaining the structure for the water molecules is not trivial. Here, the in-house code based on the genetic algorithm is used to find the global energy-minimum (see the supplementary material). The genetic-algorithm structure search is analogous to the evolutionary process in the biology. For a population consisting of a finite number of structures, the structures with lower free energies are more likely to be selected to combine into the child generation, similar to the natural selection. By iterating such selection processes from the parent-to-child generation, it is possible to find the global minimum, given enough number of generations. In this case, we add 14 water molecules around ⋆CO and hydronium. ⋆CO and hydronium are fixed during the evolutionary iterations (see the supplementary material). In Fig. 5(a), the obtained final structure is shown. To make the free energies of the initial and final structures after the reaction to be the same, four additional electrons are added to the system. Following the same procedure for the implicit solvent model case, we perform ground-state NEB to find the reaction path and reverse and tune the velocities to get a BO-MD reaction barrier. Here, the BO-MD or NEB calculated adiabatic barrier is higher than that with implicit solvents. This is because at the transition state where the proton is in the middle of CO and H₂O, unlike in the implicit solvent case, the polarization screening effect is insufficient due to the slow movement of the surrounding water molecules. rt-TDDFT is carried out with the same initial condition that is used for BO-MD. However, the NA effect becomes important near the reaction barrier similar to the implicit solvent case: rt-TDDFT and BO-MD yield opposite results for the reaction as shown in Fig. 5(b). The electronic structure’s evolution by BO-MD is illustrated in Fig. 5 of the supplementary material, including their occupations. As the reaction goes, the eigen-energy of the state is lowered, indicating the development of hybridization between the proton and ⋆CO. Meanwhile, BO-MD predicts the increased occupation of this state. But, rt-TDDFT illustrates a constant value for the occupations [shown in Fig. 5(c) of the supplementary material], although the energies of the adiabatic states is lowered owing to the hybridization. Eventually, the reaction does not happen, and it returns back to the initial structure.

In the abovementioned scenario, the number of added electrons to the whole system (mostly to the Cu layers) is used to mimic the applied overpotential. Furthermore, we perform the calculation with three additional electrons. In this case, the free energy of the final ⋆COH structure is higher than the initial hydronium state by about 0.2 eV. The calculated ground-state barrier is increased from 0.36 eV to 0.41 eV. Using BO-MD and rt-TDDFT, the reaction paths and the evolutions of the electronic structure including occupations are shown in Fig. 6 of the supplementary material. These results indicate clearly that the NA effect still plays a role when the applied potential is altered. We also perform both spin-polarized and spin-unpolarized calculations. These two types of calculation give almost the same reaction path and eigen-energy/occupation change during the reaction, i.e., the system is always non-magnetic during the reaction. Table I lists the reaction barrier calculated by the adiabatic methods and rt-TDDFT involving non-adiabaticity (using higher initial velocity to make the reaction happen). From Table I, we see that although the three cases (implicit solvent model, explicit solvent model, and different overpotentials) have rather different barriers, the barrier increase due to NA effect are all similar, around 10 meV.

To summarize, using one step of the CO₂ reduction on the copper (111) surface (⋆CO + H₃O⁺ + e⁻ → ⋆COH + H₂O) as an example, we investigated the role of the NA effect in influencing the reaction. In this reaction, the proton of hydronium influences ⋆CO to form ⋆COH. By tuning the initial velocity of the proton and monitoring the reaction using ground-state BO-MD, we identified the adiabatic reaction barrier to be the lowest initial kinetic energy of the proton that enables the reaction to happen. Using the same initial kinetic energy and structure, we observed that although BO-MD can finish the reaction, rt-TDDFT simulation with the NA effect does not induce the reaction, instead, returns the proton back to hydronium. A higher kinetic energy must be supplied to drive the proton move over the barrier to form the final structure in rt-TDDFT simulation. We also found different pictures in terms of state occupation between the BO-MD induced reaction and rt-TDDFT induced reaction. Their occupations of two CO–H related bond states are different, despite the fact that the spatial charge transfer measured from a x-y horizontal plane seem to be the same. Additional electrons are added to the system to mimic the applied overpotential to the electrode. Both implicit continuum solvent and explicit water solvent are used to simulate the same reaction. However, the NA effect remains in all the cases. Our calculation demonstrates that the NA effect increases the reaction barrier by about 10 meV for all the models that have been tested.

The supplementary material includes MD simulation under different time steps, a description of genetic algorithm, simulations with spin, and BO-MD and TDDFT simulations for explicit water cases.

This study is based on the work performed by the Joint Center for Artificial Photosynthesis, a DOE Energy Innovation Hub, supported through the Office of Science of the U.S. Department of Energy under Award No. DE-SC0004993. We use the resource of the National Energy Research Scientific Computing Center (NERSC) located in Lawrence Berkeley National Laboratory and the computational resource of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory under the Innovative and Novel Computational Impact on Theory and Experiment project.