Non-adiabatic chemical reaction refers to the electronic excitation during reactions. This effect cannot be modeled by the ground-state Born–Oppenheimer molecular dynamics (BO-MD), where the electronic structure is at the ground state for every step of ions’ movement. Although the non-adiabatic effect has been explored extensively in gas phase reactions, its role in electrochemical reactions, such as water splitting and CO2 reduction, in electrolyte has been rarely explored. On the other hand, electrochemical reactions usually involve electron transport; thus, a non-adiabatic process can naturally play a significant role. In this work, using one-step CO2 reduction as an example, we investigated the role of the non-adiabatic effect in the reaction. The reaction barriers were computed by adiabatic BO-MD and non-adiabatic real-time time dependent density functional theory (rt-TDDFT). We found that by including the non-adiabatic effect, rt-TDDFT could increase the reaction barrier up to 6% compared to the BO-MD calculated barrier when the solvent model is used to represent water. Simulations were carried out using explicit water molecules around the reaction site under different overpotentials, and similar non-adiabatic effects were found.

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 evolution1–5 and hydrogen evolution reactions,6–9 CO2 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,17 nudged elastic bands (NEBs),18 and Born–Oppenheimer molecular dynamics (BO-MD),12 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,19 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,26 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 O2 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 N2 on Ru(001), further showed the NA effect from an ab initio simulation.38 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 H2/Cu (110) and N2/W (110). However, their simulation has shown a marginal effect of non-adiabaticity on the adsorption process of diatomic molecules.37 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 CO2 adsorption to various metal surfaces has shown very fast electronic hybridization compared to the adsorption,21 exhibiting the adiabatic nature of the chemical adsorption process. Different from the initial adsorption process for CO2 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 CO2 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 CO2 reduction ⋆CO + H3O+ + e → ⋆COH + H2O (⋆: 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,42 proton-assisted chemical reactions,43 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 tn with the time step Δt = tn+1tn ≤ 0.1 fs. The time-dependent wave function is expanded as follows:

ψi(t)=jCi,j(t)ϕj(t),

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 Cij(t) for the wave function is evolved from tn to tn+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(tn) to H(tn+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 tn to tn+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 H3O+ 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 PWmat45,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 pseudopotentials47 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)50 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 CO2 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.51 Using hydronium as an example, the equilibrium reaction H3O+(aq)+e12H2(g)+H2O(aq) happens at potential U = 0 V. Thus, the enthalpy of H3O+(aq) can be expressed as H(H3O+(aq)) = 1/2H(H2(g)) + H(H2O(aq)) + 4.44 eV. Here, H stands for enthalpy, and H(H2O(aq)) − H(H2O(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 H3O+(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(H3O+(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 + H3O+ + e → ⋆COH + H2O, the energy difference ΔE = EfinalEinitial 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 H3O+ 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 H2O. 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.

FIG. 1.

Atomic structures of the (a) initial and (b) final structures. Top: top view; bottom: side view of the two structures. During the reaction, one proton of hydronium moves from O of H3O+ to O of ⋆CO. Golden: Cu; red: O; brown: C; light violet: H.

FIG. 1.

Atomic structures of the (a) initial and (b) final structures. Top: top view; bottom: side view of the two structures. During the reaction, one proton of hydronium moves from O of H3O+ to O of ⋆CO. Golden: Cu; red: O; brown: C; light violet: H.

Close modal
FIG. 2.

Reaction paths computed by adiabatic BO-MD and NA rt-TDDFT. It records the distance of the proton to oxygen of hydronium and the proton to oxygen of ⋆CO. If the reaction proceeds, the black line and red line switch, indicating the proton transferring from hydronium to ⋆CO. Otherwise, these two lines will return back.

FIG. 2.

Reaction paths computed by adiabatic BO-MD and NA rt-TDDFT. It records the distance of the proton to oxygen of hydronium and the proton to oxygen of ⋆CO. If the reaction proceeds, the black line and red line switch, indicating the proton transferring from hydronium to ⋆CO. Otherwise, these two lines will return back.

Close modal

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 H3O+. 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.

FIG. 3.

(a) Eigen energies and occupations of the states near the Fermi energy, extracted from the BO-MD simulation. (b) I: Charge density of the state at time t = 0 fs with eigen energy around −0.01 eV. Its initial occupation is 1.1. II: Charge density for the state with eigen-energy around −0.04 eV at t = 0 fs. Its initial occupation is around 1.45. (c) Eigen-energies and occupations of the states near the Fermi energy, extracted from rt-TDDFT simulation with the same initial structure and velocity to BO-MD. (d) Occupations of the adiabatic state (I and II) as a function of time for BO-MD and rt-TDDFT simulations.

FIG. 3.

(a) Eigen energies and occupations of the states near the Fermi energy, extracted from the BO-MD simulation. (b) I: Charge density of the state at time t = 0 fs with eigen energy around −0.01 eV. Its initial occupation is 1.1. II: Charge density for the state with eigen-energy around −0.04 eV at t = 0 fs. Its initial occupation is around 1.45. (c) Eigen-energies and occupations of the states near the Fermi energy, extracted from rt-TDDFT simulation with the same initial structure and velocity to BO-MD. (d) Occupations of the adiabatic state (I and II) as a function of time for BO-MD and rt-TDDFT simulations.

Close modal

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 [fj(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 fj(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 H3O+ 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.

FIG. 4.

(a) Reaction paths for BO-MD and rt-TDDFT. Both simulations have proton bonding to ⋆CO to make the reaction successful (thus rt-TDDFT has a higher initial velocity than BO-MD). (b) Measured change of total charge counted above the plane. This horizontal plane has its z-value in the middle, between C and the top Cu layer. Here, simulation- “TDDFT reaction fail” has the same initial velocity to “BO-MD reaction successful,” while “TDDFT reaction successful” has higher initial velocity than “BO-MD reaction successful.” (c) Eigen-energy of the adiabatic states for “BO-MD reaction successful” and “TDDFT reaction successful.” The bottom isosurface is the state I charge density difference between BO-MD and rt-TDDFT at t = 20 fs (charge density at the “red” dot minus the “blue” dot). In the isosurface, yellow indicates positive and blue indicates negative.

FIG. 4.

(a) Reaction paths for BO-MD and rt-TDDFT. Both simulations have proton bonding to ⋆CO to make the reaction successful (thus rt-TDDFT has a higher initial velocity than BO-MD). (b) Measured change of total charge counted above the plane. This horizontal plane has its z-value in the middle, between C and the top Cu layer. Here, simulation- “TDDFT reaction fail” has the same initial velocity to “BO-MD reaction successful,” while “TDDFT reaction successful” has higher initial velocity than “BO-MD reaction successful.” (c) Eigen-energy of the adiabatic states for “BO-MD reaction successful” and “TDDFT reaction successful.” The bottom isosurface is the state I charge density difference between BO-MD and rt-TDDFT at t = 20 fs (charge density at the “red” dot minus the “blue” dot). In the isosurface, yellow indicates positive and blue indicates negative.

Close modal

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 H2O–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)ψBOMD2(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 ψBOMD2(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.55 

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:56 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 H2O, 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.

FIG. 5.

(a) Structure of Cu/CO/H3O+ with additional 14 H2O molecules around the reaction site. Top: top view, Bottom: side view. Black and green dashed lines are hydrogen bonds between water molecules and hydrogen bonds between hydronium and ⋆CO, respectively. Hydronium ions are highlighted with different colors (violet: oxygen; orange: hydrogen) for clarity. (b) Reaction paths simulated by BO-MD and rt-TDDFT with same initial structures and velocities. Similar to Fig. 2, it records the distance of the proton to oxygen of hydronium and the proton to oxygen of ⋆CO.

FIG. 5.

(a) Structure of Cu/CO/H3O+ with additional 14 H2O molecules around the reaction site. Top: top view, Bottom: side view. Black and green dashed lines are hydrogen bonds between water molecules and hydrogen bonds between hydronium and ⋆CO, respectively. Hydronium ions are highlighted with different colors (violet: oxygen; orange: hydrogen) for clarity. (b) Reaction paths simulated by BO-MD and rt-TDDFT with same initial structures and velocities. Similar to Fig. 2, it records the distance of the proton to oxygen of hydronium and the proton to oxygen of ⋆CO.

Close modal

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.

TABLE I.

Reaction barriers calculated by NEB, BO-MD, and rt-TDDFT. Here, NEB and BO-MD are only ground-state calculations, whereas rt-TDDFT involves the NA effect beyond the ground-state approximation. The last column is the percentage change of the barrier by the NA effect.

Reaction barrier EENEB (eV)EBO-MD (eV)Ert-TDDFT (eV)Ert-TDDFTEBO-MD (meV)
Implicit solvent (add 2e0.080 0.196 0.208 12 
Hybrid solvent (add 4e0.360 0.288 0.299 11 
Hybrid solvent (add 3e0.411 0.363 0.373 10 
Reaction barrier EENEB (eV)EBO-MD (eV)Ert-TDDFT (eV)Ert-TDDFTEBO-MD (meV)
Implicit solvent (add 2e0.080 0.196 0.208 12 
Hybrid solvent (add 4e0.360 0.288 0.299 11 
Hybrid solvent (add 3e0.411 0.363 0.373 10 

To summarize, using one step of the CO2 reduction on the copper (111) surface (⋆CO + H3O+ + e → ⋆COH + H2O) 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.

1.
Á.
Valdés
,
Z.-W.
Qu
,
G.-J.
Kroes
,
J.
Rossmeisl
, and
J. K.
Nørskov
, “
Oxidation and photo-oxidation of water on TiO2 surface
,”
J. Phys. Chem. C
112
,
9872
9879
(
2008
).
2.
T. A.
Pham
,
Y.
Ping
, and
G.
Galli
, “
Modelling heterogeneous interfaces for solar water splitting
,”
Nat. Mater.
16
,
401
408
(
2017
).
3.
H. H.
Pham
,
M.-J.
Cheng
,
H.
Frei
, and
L.-W.
Wang
, “
Surface proton hopping and fast-kinetics pathway of water oxidation on Co3O4(001) surface
,”
ACS Catal.
6
,
5610
5617
(
2016
).
4.
Y.
Wu
,
P.
Lazic
,
G.
Hautier
,
K.
Persson
, and
G.
Ceder
, “
First principles high throughput screening of oxynitrides for water-splitting photocatalysts
,”
Energy Environ. Sci.
6
,
157
168
(
2013
).
5.
Y.
Zhou
,
G.
Gao
,
Y.
Li
,
W.
Chu
, and
L.-W.
Wang
, “
Transition-metal single atoms in nitrogen-doped graphenes as efficient active centers for water splitting: A theoretical study
,”
Phys. Chem. Chem. Phys.
21
,
3024
3032
(
2019
).
6.
J.
Greeley
,
T. F.
Jaramillo
,
J.
Bonde
,
I.
Chorkendorff
, and
J. K.
Nørskov
, “
Computational high-throughput screening of electrocatalytic materials for hydrogen evolution
,”
Nat. Mater.
5
,
909
913
(
2006
).
7.
M.
Chhetri
,
S.
Maitra
,
H.
Chakraborty
,
U. V.
Waghmare
, and
C. N. R.
Rao
, “
Superior performance of borocarbonitrides, BxCyNz, as stable, low-cost metal-free electrocatalysts for the hydrogen evolution reaction
,”
Energy Environ. Sci.
9
,
95
101
(
2016
).
8.
Q.
Tang
and
D.-E.
Jiang
, “
Mechanism of hydrogen evolution reaction on 1T-MoS2 from first principles
,”
ACS Catal.
6
,
4953
4961
(
2016
).
9.
G.
Gao
,
A. P.
O’Mullane
, and
A.
Du
, “
2D MXenes: A new family of promising catalysts for the hydrogen evolution reaction
,”
ACS Catal.
7
,
494
500
(
2017
).
10.
A. A.
Peterson
,
F.
Abild-Pedersen
,
F.
Studt
,
J.
Rossmeisl
, and
J. K.
Nørskov
, “
How copper catalyzes the electroreduction of carbon dioxide into hydrocarbon fuels
,”
Energy Environ. Sci.
3
,
1311
1315
(
2010
).
11.
J. H.
Montoya
,
A. A.
Peterson
, and
J. K.
Nørskov
, “
Insights into C–C coupling in CO2 electroreduction on copper electrodes
,”
ChemCatChem
5
,
737
742
(
2013
).
12.
T.
Cheng
,
H.
Xiao
, and
W. A.
Goddard
, “
Free-energy barriers and reaction mechanisms for the electrochemical reduction of CO on the Cu(100) surface, including multiple layers of explicit solvent at pH 0
,”
J. Phys. Chem. Lett.
6
,
4767
4773
(
2015
).
13.
T.
Cheng
,
H.
Xiao
, and
W. A.
Goddard
, “
Reaction mechanisms for the electrochemical reduction of CO2 to CO and formate on the Cu(100) surface at 298 K from quantum mechanics free energy calculations with explicit water
,”
J. Am. Chem. Soc.
138
,
13802
13805
(
2016
).
14.
Y.
Ping
,
W. A.
Goddard
, and
G. A.
Galli
, “
Energetics and solvation effects at the photoanode/catalyst interface: Ohmic contact versus Schottky barrier
,”
J. Am. Chem. Soc.
137
,
5264
5267
(
2015
).
15.
A. G.
Scheuermann
,
J. P.
Lawrence
,
K. W.
Kemp
,
T.
Ito
,
A.
Walsh
,
C. E. D.
Chidsey
,
P. K.
Hurley
, and
P. C.
McIntyre
, “
Design principles for maximizing photovoltage in metal-oxide-protected water-splitting photoanodes
,”
Nat. Mater.
15
,
99
105
(
2016
).
16.
Q.
Yan
,
J.
Yu
,
S. K.
Suram
,
L.
Zhou
,
A.
Shinde
,
P. F.
Newhouse
,
W.
Chen
,
G.
Li
,
K. A.
Persson
,
J. M.
Gregoire
, and
J. B.
Neaton
, “
Solar fuels photoanode materials discovery by integrating high-throughput theory and experiment
,”
Proc. Natl. Acad. Sci. U. S. A.
114
,
3040
3043
(
2017
).
17.
J. K.
Nørskov
,
J.
Rossmeisl
,
A.
Logadottir
,
L.
Lindqvist
,
J. R.
Kitchin
,
T.
Bligaard
, and
H.
Jónsson
, “
Origin of the overpotential for oxygen reduction at a fuel-cell cathode
,”
J. Phys. Chem. B
108
,
17886
17892
(
2004
).
18.
B. J.
Berne
,
G.
Ciccotti
, and
D. F.
Coker
,
Classical and Quantum Dynamics in Condensed Phase Simulations
(
Villa Marigola
,
Lerici
,
1998
).
19.
J.
Rossmeisl
,
K.
Chan
,
R.
Ahmed
,
V.
Tripković
, and
M. E.
Björketun
, “
pH in atomic scale simulations of electrochemical interfaces
,”
Phys. Chem. Chem. Phys.
15
,
10321
10325
(
2013
).
20.
S.
Ringe
,
E. L.
Clark
,
J.
Resasco
,
A.
Walton
,
B.
Seger
,
A. T.
Bell
, and
K.
Chan
, “
Understanding cation effects in electrochemical CO2 reduction
,”
Energy Environ. Sci.
12
,
3001
3014
(
2019
).
21.
M.
Bajdich
,
M.
Fields
,
L. D.
Chen
,
R. B.
Sandberg
,
K.
Chan
, and
J. K.
Nørskov
, “
Electron transfer to CO2 during adsorption at the metal–solution interface
,”
Phys. Chem. C
123
,
29278
(
2019
).
22.
O.
Andreussi
,
I.
Dabo
, and
N.
Marzari
, “
Revised self-consistent continuum solvation in electronic-structure calculations
,”
J. Chem. Phys.
136
,
064102
(
2012
).
23.
K.
Mathew
,
R.
Sundararaman
,
K.
Letchworth-Weaver
,
T. A.
Arias
, and
R. G.
Hennig
, “
Implicit solvation model for density-functional study of nanocrystal surfaces and reaction pathways
,”
J. Chem. Phys.
140
,
084106
(
2014
).
24.
G. A.
Tritsaris
,
D.
Vinichenko
,
G.
Kolesov
,
C. M.
Friend
, and
E.
Kaxiras
, “
Dynamics of the photogenerated hole at the rutile TiO2(110)/water interface: A nonadiabatic simulation study
,”
J. Phys. Chem. C
118
,
27393
27401
(
2014
).
25.
S.
Hammes-Schiffer
, “
Theory of proton-coupled electron transfer in energy conversion processes
,”
Acc. Chem. Res.
42
,
1881
1889
(
2009
).
26.
Z.
Wang
,
S.-S.
Li
, and
L.-W.
Wang
, “
Efficient real-time time-dependent density functional theory method and its application to a collision of an ion with a 2D material
,”
Phys. Rev. Lett.
114
,
063004
(
2015
).
27.
G.
Bi
,
J.
Kang
, and
L.-W.
Wang
, “
High velocity proton collision with liquid lithium: A time dependent density functional theory study
,”
Phys. Chem. Chem. Phys.
19
,
9053
9058
(
2017
).
28.
D. C.
Yost
,
Y.
Yao
, and
Y.
Kanai
, “
Examining real-time time-dependent density functional theory nonequilibrium simulations for the calculation of electronic stopping power
,”
Phys. Rev. B
96
,
115134
(
2017
).
29.
D.
Wang
,
Z.-P.
Liu
, and
W.-M.
Yang
, “
Proton-promoted electron transfer in photocatalysis: Key step for photocatalytic hydrogen evolution on metal/titania composites
,”
ACS Catal.
7
,
2744
2752
(
2017
).
30.
M.
Alducin
,
R.
Díez Muiño
, and
J. I.
Juaristi
, “
Non-adiabatic effects in elementary reaction processes at metal surfaces
,”
Prog. Surf. Sci.
92
,
317
340
(
2017
).
31.
G.-J.
Kroes
,
J. I.
Juaristi
, and
M.
Alducin
, “
Vibrational excitation of H2 scattering from Cu(111): Effects of surface temperature and of allowing energy exchange with the surface
,”
J. Phys. Chem. C
121
,
13617
13633
(
2017
).
32.
B.
Jiang
,
M.
Alducin
, and
H.
Guo
, “
Electron–hole pair effects in polyatomic dissociative chemisorption: Water on Ni(111)
,”
J. Phys. Chem. Lett.
7
,
327
331
(
2016
).
33.
X.
Luo
,
B.
Jiang
,
J. I.
Juaristi
,
M.
Alducin
, and
H.
Guo
, “
Electron-hole pair effects in methane dissociative chemisorption on Ni(111)
,”
J. Chem. Phys.
145
,
044704
(
2016
).
34.
G.
Füchsel
,
M.
del Cueto
,
C.
Díaz
, and
G.-J.
Kroes
, “
Enigmatic HCl + Au(111) reaction: A puzzle for theory and experiment
,”
J. Phys. Chem. C
120
,
25760
25779
(
2016
).
35.
I.
Goikoetxea
,
J. I.
Juaristi
,
M.
Alducin
, and
R.
Díez Muiño
, “
Dissipative effects in the dynamics of N2 on tungsten surfaces
,”
J. Phys.: Condens. Matter
21
,
264007
(
2009
).
36.
C.
Carbogno
,
J.
Behler
,
K.
Reuter
, and
A.
Groß
, “
Signatures of nonadiabatic O2 dissociation at Al(111): First-principles fewest-switches study
,”
Phys. Rev. B
81
,
035410
(
2010
).
37.
J. I.
Juaristi
,
M.
Alducin
,
R. D.
Muiño
,
H. F.
Busnengo
, and
A.
Salin
, “
Role of electron-hole pair excitations in the dissociative adsorption of diatomic molecules on metal surfaces
,”
Phys. Rev. Lett.
100
,
116102
(
2008
).
38.
L.
Diekhöner
,
L.
Hornekær
,
H.
Mortensen
,
E.
Jensen
,
A.
Baurichter
,
V. V.
Petrunin
, and
A. C.
Luntz
, “
Indirect evidence for strong nonadiabatic coupling in N2 associative desorption from and dissociative adsorption on Ru(0001)
,”
J. Chem. Phys.
117
,
5018
5030
(
2002
).
39.
S. A.
Fischer
,
W. R.
Duncan
, and
O. V.
Prezhdo
, “
Ab initio nonadiabatic molecular dynamics of wet-electrons on the TiO2 surface
,”
J. Am. Chem. Soc.
131
,
15483
15491
(
2009
).
40.
H.
Petek
and
J.
Zhao
, “
Ultrafast interfacial proton-coupled electron transfer
,”
Chem. Rev.
110
,
7082
7099
(
2010
).
41.
B. G.
Oscar
,
W.
Liu
,
N. D.
Rozanov
, and
C.
Fang
, “
Ultrafast intermolecular proton transfer to a proton scavenger in an organic solvent
,”
Phys. Chem. Chem. Phys.
18
,
26151
26160
(
2016
).
42.
N. J.
Turro
,
Modern Molecular Photochemistry
(
University Science Books
,
1991
).
43.
S.
Hammes-Schiffer
, “
Theoretical perspectives on proton-coupled electron transfer reactions
,”
Acc. Chem. Res.
34
,
273
281
(
2001
).
44.
Y.
Miyamoto
and
H.
Zhang
, “
Electronic excitation in an Ar7+ ion traversing a graphene sheet: Molecular dynamics simulations
,”
Phys. Rev. B
77
,
161402
(
2008
).
45.
W.
Jia
,
Z.
Cao
,
L.
Wang
,
J.
Fu
,
X.
Chi
,
W.
Gao
, and
L.-W.
Wang
, “
The analysis of a plane wave pseudopotential density functional theory code on a GPU machine
,”
Comput. Phys. Commun.
184
,
9
18
(
2013
).
46.
W.
Jia
,
J.
Fu
,
Z.
Cao
,
L.
Wang
,
X.
Chi
,
W.
Gao
, and
L.-W.
Wang
, “
Fast plane wave density functional theory molecular dynamics calculations on multi-GPU machines
,”
J. Comput. Phys.
251
,
102
115
(
2013
).
47.
D. R.
Hamann
, “
Optimized norm-conserving Vanderbilt pseudopotentials
,”
Phys. Rev. B
88
,
085117
(
2013
).
48.
G.
Kresse
,
A.
Gil
, and
P.
Sautet
, “
Significance of single-electron energies for the description of CO on Pt(111)
,”
Phys. Rev. B
68
,
073401
(
2003
).
49.
S. E.
Mason
,
I.
Grinberg
, and
A. M.
Rappe
, “
First-principles extrapolation method for accurate CO adsorption energies on metal surfaces
,”
Phys. Rev. B
69
,
161401
(
2004
).
50.
B.
Hammer
,
L. B.
Hansen
, and
J. K.
Nørskov
, “
Improved adsorption energetics within density-functional theory using revised Perdew-Burke-Ernzerhof functionals
,”
Phys. Rev. B
59
,
7413
7421
(
1999
).
51.
V. S.
Bryantsev
,
M. S.
Diallo
, and
W. A.
Goddard
 III
, “
Calculation of solvation free energies of charged solutes using mixed cluster/continuum models
,”
J. Phys. Chem. B
112
,
9709
9719
(
2008
).
52.
C.-G.
Zhan
and
D. A.
Dixon
, “
First-principles determination of the absolute hydration free energy of the hydroxide ion
,”
J. Phys. Chem. A
106
,
9737
9744
(
2002
).
53.
C.
Dupont
,
O.
Andreussi
, and
N.
Marzari
, “
Self-consistent continuum solvation (SCCS): The case of charged systems
,”
J. Chem. Phys.
139
,
214110
(
2013
).
54.
M. W.
Palascak
and
G. C.
Shields
, “
Accurate experimental values for the free energies of hydration of H+, OH, and H3O+
,”
J. Phys. Chem. A
108
,
3692
3694
(
2004
).
55.
Y.
Quan
,
V.
Pardo
, and
W. E.
Pickett
, “
Formal valence,3d-electron occupation, and charge-order transitions
,”
Phys. Rev. Lett.
109
,
216401
(
2012
).
56.
J. A.
Gauthier
,
S.
Ringe
,
C. F.
Dickens
,
A. J.
Garza
,
A. T.
Bell
,
M.
Head-Gordon
,
J. K.
Nørskov
, and
K.
Chan
, “
Challenges in modeling electrochemical reaction energetics with polarizable continuum models
,”
ACS Catal.
9
,
920
931
(
2019
).

Supplementary data