The catalytic reduction in carbon dioxide is a crucial step in many chemical industrial reactions, such as methanol synthesis, the reverse water–gas shift reaction, and formic acid synthesis. Here, we investigate the role of bulk hydrogen, where hydrogen atoms are found deep inside a metal surface as opposed to subsurface ones, upon CO2 reduction over a Ni(110) surface using density functional theory and ab initio molecular dynamics simulations. While it has previously been shown that subsurface hydrogen stabilizes CO2 and can aid in overcoming reaction barriers, the role of bulk hydrogen is less studied and thus unknown with regard to CO2 reduction. We find that the presence of bulk hydrogen can significantly alter the electronic structure of the Ni(110) surface, particularly the work function and d-band center, such that CO2 adsorbs more strongly to the surface and is more easily reduced. Our results show an enhanced CO2 dissociation in the presence of bulk hydrogen, shedding light on a hitherto underappreciated mechanistic pathway for CO2 reduction on metal surfaces.

The catalytic reduction in carbon dioxide is a crucial step in many chemical reactions, such as methanol synthesis, the reverse water–gas shift (rWGS) reaction, the Bosch reaction, the Sabatier reaction, and formic acid synthesis.1–18 These reactions are typically catalyzed by transition metals, including Cu, Co, and Ni, of which the surface electronic structure primarily determines the overall activation energy of the reaction.19–21 Depending on reaction conditions (temperature and pressure), hydrogen will be present on the surface (Hsurf), below the surface (Hsub), and/or in the bulk (Hbulk). Under low-pressure (e.g., UHV) conditions, Hsub and Hbulk can be formed through two mechanisms: (i) individual H atoms hitting the surface with sufficient energy to overcome the energy barrier [∼1 eV is needed for a Ni(111) surface] and (ii) collision of an inert gas atom or molecule with Hsurf with a sufficient kinetic energy (∼2.6 eV) to push it into the subsurface.22,23 However, heterogeneous catalysis of hydrogenation reactions typically occurs at high pressures (typically ∼20 bar but even up to ∼100 bar),14,24 at which hydrogen more easily populates the bulk of the catalyst. The precise location (e.g., octahedral/tetrahedral site; see Figs. 1(a) and 1(b)), number of layers below surface, proximity to adsorbed species, and concentration of each type of hydrogen determine the electronic structure of the catalyst surface and can thereby substantially impact its catalytic performance. While the correlation between the Hsub content and hydrogenation activity has been investigated since the 1950s,1,25 its precise role, either as a spectator or as a reactant, is often unclear.26–30 Generally, however, it is well-known that Hsub is less stable and becomes more reactive than Hsurf when it emerges to the surface, although its reactivity depends on the barrier to emerge.2,31 When such an emergence occurs, the extra energy released helps overcome the barrier to hydrogenation of adsorbed species.32 Furthermore, Hsub tends to stabilize adsorbates by inducing surface strain that shifts the d-band of the metallic surface upward.21,33–35 Less is known regarding the role of Hbulk, which is often assumed to play a negligible role, if any at all, in catalytic hydrogenation.

FIG. 1.

Nickel face-centered cubic (fcc) crystal structure with hydrogen in (a) octahedral sites and (b) tetrahedral sites. (c) Top view and (d) side view of the Ni(110) slab used in all simulations. A primitive unit cell is shown in (c) by the red dashed lines.

FIG. 1.

Nickel face-centered cubic (fcc) crystal structure with hydrogen in (a) octahedral sites and (b) tetrahedral sites. (c) Top view and (d) side view of the Ni(110) slab used in all simulations. A primitive unit cell is shown in (c) by the red dashed lines.

Close modal

Nickel, in particular, exhibits promising catalytic activity for reactions involving the dissociation of CO2 to CO, such as the rWGS (CO2 + H2 → CO + H2O) and the dry reforming reaction (CO2 + CH4 → 2H2 + 2CO), as well as CO methanation (CO + 3H2 → CH4 + H2O) and methanol production (CO + 2H2 → CH3OH). On well-studied Ni surfaces, chemisorption of CO2 is most favorable and exothermic on Ni(110), with Ni(100) being less favorable and Ni(111) being the least favorable.36,37 Ni(110) has also been shown to be more reactive than Ni(100) and Ni(111),38 as it has a greater resemblance to reactive, coordinatively unsaturated sites. Interestingly, recent experimental and theoretical studies have demonstrated that the presence of Hsub atoms alters the electronic structure of the Ni surface in a way that significantly changes the reaction mechanism.39–41 For example, Schatz et al., using density functional theory (DFT) calculations, showed that in the case of the rWGS, the redox pathway dominates in the presence of Hsub, while the associative pathway dominates in their absence.21 They also demonstrated that Hsub atoms have the potential to aid greatly in the catalysis of CO2 dissociation by favoring dissociation over desorption. Furthermore, in the case of CO2 hydrogenation to formic acid, studies have shown that the emergence of Hsub to the surface changes the reaction from endothermic to exothermic.32,42

In general, these studies assume the presence of co-adsorbed hydrogen (i.e., Hsurf atoms) and only regard H atoms in the octahedral (Oh) sites below the first Ni layer [see Fig. 1(a)] as Hsub, claiming that other Hbulk sites do not significantly affect surface Ni atoms and adsorbed species.

In addition, the theoretical studies mentioned above utilized static DFT calculations, which do not intrinsically take into account thermal effects, although the latter can be added as a correction through the harmonic approximation. Molecular dynamics also allow us to better sample the configurational space available to the system at elevated temperatures. As the distribution of these H species is strongly influenced by the temperature, high temperatures (e.g., ∼300–400 °C in WGS) can greatly increase the diffusion processes of H atoms from the surface to the bulk and vice versa. To address these gaps in the literature, here, we decided to study the effects of Hbulk, Hsub, and Hsurf on the adsorption and reduction in CO2 on the Ni(110) surface, which is considered a representative surface model for the activity of low-coordination sites of Ni catalysts. Furthermore, using ab initio molecular dynamics (AIMD) simulations, we are able to follow how the coordination and location of the H atoms, not only on the surface but also on the sub-surface and bulk, evolve over time and how this distribution changes with temperature. While our study confirms the literature’s findings that Hsub favors CO2 dissociation and Hsurf favors CO2 desorption, we also observe that (i) Hbulk significantly alters both the structure and electronic states of the Ni surface and (ii) CO2 exothermically dissociates in the presence of only Hbulk, i.e., even in the absence of any Hsurf and Hsub. Our results indicate that there is a much more complex interplay among surface, subsurface, and bulk H than previously suspected and also allude to the possibility of other reaction pathways not discussed in previous studies.

To simulate the metal substrate, we adopted a simulation cell consisting of a 3 × 3 Ni (110) slab with seven Ni layers and 12 Å vacuum layer; see Figs. 1(c) and 1(d). Following the convention of Schatz et al., we refer to any H atoms below the second Ni layer as Hbulk and H atoms below the first and above the second Ni layer as Hsub. H atoms on the surface (above the first Ni layer) are referred to as Hsurf. Examples of Hsurf/Ni(110), Hsub/Ni(110), and Hbulk/Ni(110) are shown in Figs. 2(a)2(c). Each of these systems at a range of hydrogen coverages was used for static density functional theory (DFT) electronic structure calculations. For the AIMD simulations, systems including six H atoms (2/3 monolayer, ML) were used to create the Hbulk/Ni(110) and Hsurf/Ni(110) models to interrogate the properties of Hbulk. During all simulations, the bottom Ni layer was fixed at the bulk positions. To produce initial structures with only Hsurf, three H2 molecules were initially placed randomly on the Ni surface and were subsequently relaxed. To produce initial structures with Hbulk and Hsub, H atoms were placed in the tetrahedral (Td) sites below the second Ni layer and the octahedral (Oh) sites below the first Ni layer, respectively; see Figs. 2(b) and 2(c). For studying the effect of H on adsorbed CO2, a single CO2 molecule was placed in the well-known chemisorbed bent CO2δ state on the “short bridge–long bridge” site, which has been shown to be the most stable adsorption site for CO2 on Ni(110).21 

FIG. 2.

Side view of (a) Hsurf/Ni(110), (b) Hsub/Ni(110), (c) Hbulk/Ni(110), (d) Hsurf/Ni(110) with CO2, (e) Hsub/Ni(110) with CO2, and (f) Hbulk/Ni(110) with CO2. Color codes: Ni in blue and H in dark white.

FIG. 2.

Side view of (a) Hsurf/Ni(110), (b) Hsub/Ni(110), (c) Hbulk/Ni(110), (d) Hsurf/Ni(110) with CO2, (e) Hsub/Ni(110) with CO2, and (f) Hbulk/Ni(110) with CO2. Color codes: Ni in blue and H in dark white.

Close modal

AIMD simulations were carried out using the CP2K program.43 The potential energy was calculated at the PBE-D3 level.44,45 The nuclei and core electrons were described through the norm-conserving Goedecker–Teter–Hutter (GTH) pseudopotentials.46,47 We adopted the Gaussian-plane wave hybrid approach48 in which the Gaussian basis sets of double-zeta valence with polarization (DZVP)49 quality in conjunction with a plane wave cut-off 400 Ry were used. Spin polarization was applied for all AIMD and static DFT calculations.

The AIMD simulations were carried out using the NVT ensemble50 for different temperatures at 100, 200, 400, and 1000 K maintained by the Nose–Hoover thermostat.51,52 Each system was first optimized to a local minimum, followed by a 10–15 ps equilibration run at a given temperature. Finally, a production run of approximately 20 ps was used for the analysis of the distribution of H atoms and other properties. All electronic structure properties (Bader charges,53–57 work functions, and projected density of states) were computed from trajectories at a coverage of 6H and from randomized configuration at other coverages. For the electronic structure calculations, we utilized a larger basis set of triple zeta quality for a more accurate description of electronic density of states and work function.

We first examine the effects of H coverage on the surface structure of Ni(110) and the electronic structure of the various H/Ni(110) systems. The Ni–Ni bond length of the surface atoms, rNi–Ni, and the Ni–Ni coordination number of the surface atoms, CNNi–Ni, of the optimized pristine Ni(110) are 2.45 Å and 7.0, respectively. rNi–Ni is within a 1% error of the experimental crystal structure Ni–Ni distance of 2.48 Å,58 while CNNi–Ni is the expected coordination number for an fcc (110) surface. In Table I, we show rNi–Ni and CNNi–Ni for the optimized Hsurf/Ni(110), Hsub/Ni(110), and Hbulk/Ni(110) systems as a function of H coverage. At low coverages, rNi–Ni is essentially unaffected, while at a coverage of 6H (2/3 ML), the bond length has increased for all H-types, with the largest increase for 6Hsub/Ni(110) of 0.07 Å relative to pristine Ni(110). We see that CNNi–Ni is also largely unaffected, with slight increases at every coverage, due to some surface Ni atoms becoming over-coordinated. In the last column of Table I, we have also tabulated the interlayer relaxation between the surface layer and the subsurface layer, Δd12, with respect to the interlayer spacing in pristine Ni(110), defined as the difference in the average height between the atoms in the first and second layers of the Ni slab. We observe that all H-types cause layer relaxation (predominantly expansion), even at low coverages. However, with increasing coverage, Hsub produces the largest change in interlayer spacing, with a change of 17.7% for 6Hsub/Ni(110), while 6Hsurf/Ni(110) and 6Hbulk/Ni(110) exhibit only a ∼6% change.

TABLE I.

Ni–Ni bond lengths, rNi–Ni, and coordination numbers, CNNi–Ni, of surface Ni atoms and layer relaxation Δd12 with respect to the pure Ni(110) interlayer spacing of the optimized structures.

rNi–Ni (Å)CNNi–NiΔd12 (%)
Ni(110) 2.45 7.0  
2Hsurf/Ni(110) 2.48 7.4 3.50 
2Hsub/Ni(110) 2.47 7.0 2.01 
2Hbulk/Ni(110) 2.47 7.2 4.25 
4Hsurf/Ni(110) 2.45 7.2 −1.5 
4Hsub/Ni(110) 2.49 7.2 9.6 
4Hbulk/Ni(110) 2.44 7.2 −5.22 
6Hsurf/Ni(110) 2.48 7.0 6.5 
6Hsub/Ni(110) 2.52 7.0 17.7 
6Hbulk/Ni(110) 2.50 7.2 6.11 
rNi–Ni (Å)CNNi–NiΔd12 (%)
Ni(110) 2.45 7.0  
2Hsurf/Ni(110) 2.48 7.4 3.50 
2Hsub/Ni(110) 2.47 7.0 2.01 
2Hbulk/Ni(110) 2.47 7.2 4.25 
4Hsurf/Ni(110) 2.45 7.2 −1.5 
4Hsub/Ni(110) 2.49 7.2 9.6 
4Hbulk/Ni(110) 2.44 7.2 −5.22 
6Hsurf/Ni(110) 2.48 7.0 6.5 
6Hsub/Ni(110) 2.52 7.0 17.7 
6Hbulk/Ni(110) 2.50 7.2 6.11 

Next, we examine how the H-types and coverage affect the surface electronic structure of Ni(110). A Bader analysis53–57 on the hydrogen and surface Ni atoms revealed that the charge per hydrogen atom, qH, does not vary significantly with respect to both type and coverage, as shown in Fig. 3(a) (the range of the y axis is only 0.04e). However, the H-type and coverage significantly alter the Ni atom charge in the first layer, qNi, as shown in Fig. 3(b) for Hsurf/Ni(110) and Hsub/Ni(110). As we would expect, qNi increases with increasing H coverage for all H-types, but most drastically for Hsurf, as plotted in Fig. 3(b). For Hbulk/Ni(110), qNi remains close to the charge of pure Ni(110), consistent with previous findings by Schatz et al.,21 with only a slight increase with coverage. In Figs. 3(c) and 3(d), we also show the work function ϕ of each system, computed according to ϕ = EvacEF, where EF is the Fermi level and Evac is the vacuum potential energy, as well as the d-band center ϵd, given by

ϵd=ρdEEdEρdEdE,

where ρd(E) is the density of states projected onto the Ni d-orbitals of the surface atoms [see Fig. S2 for ρd(E) of the selected systems]. Each describes a different catalytic property of the surface, where ϕ is an indicator of the reducibility59 of the surface and ϵd is the stability of adsorbates.60 We see that, generally, all H-types reduce the work function of Ni(110) as most values are below the black dashed line in Fig. 3(c).

FIG. 3.

(a) Charge per hydrogen atom qH, (b) total surface charge qNi, (c) work function ϕ, and (d) d-band center ϵd of H/Ni(110) systems as a function of H coverage. The black dashed lines represent the corresponding value of pristine Ni(110).

FIG. 3.

(a) Charge per hydrogen atom qH, (b) total surface charge qNi, (c) work function ϕ, and (d) d-band center ϵd of H/Ni(110) systems as a function of H coverage. The black dashed lines represent the corresponding value of pristine Ni(110).

Close modal

Interestingly, Hsub/Ni(110) exhibits an increase in ϕ with coverage, while Hbulk/Ni(110) shows a decrease with coverage. Of all Hbulk systems under study, 6Hbulk/Ni(110) has the lowest work function (4.49 eV) and is expected to be the strongest reducing agent. Now, examining Fig. 3(d), we see that the d-band center becomes less negative with a coverage as beyond 0.5 ML. Typically, such an increase in ϵd indicates increased adsorbate stability,35 which has been well-documented in the literature for Hsub/Ni(110).33 While our results also show that ϵd is coverage dependent, they imply that higher coverages increase adsorbate stability. On the other hand, ϕ for all H/Ni(110) systems is lower than the value of pure Ni(110), which is desirable for reactions involving CO2 reduction (e.g., rWGS). We will explore the effects of the interplay of ϵd and ϕ in these systems later in the text.

Next, we examine hydrogen diffusion and thermal fluctuations of Ni coordination at a range of temperatures from AIMD simulations for 6Hsurf/Ni(110) and 6Hbulk/Ni(110). We chose these systems because of the lack of discussion on Hbulk in the literature, and we chose a coverage of 2/3 ML to ensure reactivity of the surface. Furthermore, we have demonstrated that at high coverages, Hbulk is promising for improving reducibility (low ϕ) and adsorbate stability (high ϵd) of Ni(110). The diffusion coefficient D can be computed from the mean squared displacement (MSD), given by

MSD=ΔR2τH=1N0NHj=1N0i=1NHRitj+τRitj2,

where Ri is the position of atom i, NH is the number of H atoms, and N0 is the number of time origins, spaced by τ. In Fig. S3, we plot the MSD along the plane, MSDxy, and normal to the surface, MSDz. We observe that Hsurf starts diffusing at 200 K, while Hbulk remains in the bulk until 400 K, where diffusion starts but still stays lower than Hsurf. We obtain the diffusion coefficients Dxy and Dz according to MSD = 4 and MSD = 2, respectively, over the linear portions of each curve [(τ = 2.5 − 6 ps); see Fig. S3], plotted in Fig. 4. Not surprisingly, we see that Dxy for Hbulk is significantly lower, by almost two orders of magnitude, than Dxy for Hsurf. Dz is much lower than Dxy for both Hbulk and Hsurf, although Hbulk is slightly more mobile than Hsurf along z.

FIG. 4.

Log of the diffusion coefficients D of Hsurf and Hbulk parallel (Dxy) and perpendicular (Dz) to the Ni surface plane as a function of 1/T, for T = 100, 200, and 400 K, with ideal Arrhenius lines. The computed activation energies are as follows: (Dxy)Hsurf: 2.12 kcal/mol, (Dz)Hsurf: 0.67 kcal/mol, (Dxy)Hbulk: 0.90 kcal/mol, and (Dz)Hbulk: 0.31 kcal/mol.

FIG. 4.

Log of the diffusion coefficients D of Hsurf and Hbulk parallel (Dxy) and perpendicular (Dz) to the Ni surface plane as a function of 1/T, for T = 100, 200, and 400 K, with ideal Arrhenius lines. The computed activation energies are as follows: (Dxy)Hsurf: 2.12 kcal/mol, (Dz)Hsurf: 0.67 kcal/mol, (Dxy)Hbulk: 0.90 kcal/mol, and (Dz)Hbulk: 0.31 kcal/mol.

Close modal

To understand the effect of temperature on the surface structure of H/Ni(110), in Figs. 5(a) and 5(b), we plot the Ni–Ni coordination number, CNNi–Ni, of each surface Ni atom over time. At 200 K and below, CNNi–Ni of all Ni surface atoms is the same as those of the pure Ni(110) surface atoms, CNNi–Ni = 7. At 400 K, we start to observe over-coordinated Ni atoms in both systems, predominantly over-coordinated (CNNi–Ni > 7) sites, especially for Hsurf/Ni(110). Simultaneously, we also see a change in the coordination of Hsurf to the surface with temperature, with a shift from hollow (threefold coordinated sites) to bridge (twofold coordinated sites), as shown in Fig. 5(c). As shown in Fig. 5(d), a Hsurf atom goes as far as pulling adjacent Ni atoms outward on the bridge site, thus causing an interlayer expansion between the top two layers of Ni(110). Through inspection of the trajectory, we also observe that the over-coordinated Ni atoms shown in Figs. 5(a) and 5(b) are due to bond formations between Ni atoms in the surface layer (layer 1). Such changes in surface geometry are expected for Hsurf/Ni(110) as Ni(110) is the most reactive of the Ni low-index surfaces. As for Hbulk/Ni(110), we observe fewer over-coordinated Ni atoms on the surface, but we note that at 400 K, one Hbulk atom diffuses to the surface within 20 picoseconds of simulation (see Fig. S4a). We attribute the changes in Ni coordination to the presence of this single Hsurf atom in the Hbulk/Ni(110) 400 K simulation. From this detail, we can obtain a crude estimate of the kinetic barrier for Hbulk to emerge to the surface, computed from the following equation:61 

k=kbTheEakbT(HbulkHsurf),

where k is the rate constant, Ea is the activation energy, kb is the Boltzmann constant, h is the Planck constant, and T is the temperature. Here, we computed Ea for k = 1 and k = 1/100 ps−1 at T = 400 K since we observed a single event within the first picosecond for a simulation time of up to ∼15 ps. From this equation, we obtain a range for the Hbulk → Hsurf free energy barrier of 0.1–0.3 eV for Hbulk to emerge to the surface, which is comparable to the reported energy barrier for Hsub to emerge, which is on the order of 0.11–0.14 eV, as estimated from static calculations.21,31 To check on the magnitude of these values, we analyzed a trajectory at a higher temperature (1000 K) at which all Hbulk diffuse to the surface during the simulation. From this simulation, we obtain a value (∼0.3 eV) comparable in magnitude to the expected range. Due to the instability of the tetrahedral Hbulk sites relative to the octahedral Hsub sites, we anticipate that the energy gained moving first from the bulk to the subsurface enables Hbulk to overcome the barrier to move from the subsurface to the surface. However, in order to fully exploit the enhanced properties of Hbulk/Ni(110), it is important to keep H in the bulk and minimize the amount of Hsurf and Hsub. Development of a route to producing Hbulk, which usually only occurs once the surface is saturated,62 could be transformative for the catalysis community.

FIG. 5.

Ni–Ni coordination number CNNi–Ni as a function of time for (a) 6Hbulk/Ni(110) and (b) 6Hsurf/Ni(110). Each curve of the same temperature corresponds to one Ni atom, i.e., the number of markers at each time t is the total number of Ni nearest neighbors. (c) Distribution of Hsurf adsorption sites at 100 K, 200 K, and 400 K for 6Hsurf/Ni(110). (d) Examples of bridge (twofold coordinated) and hollow (threefold coordinated) sites for Hsurf shown with the orange and red bonds, respectively.

FIG. 5.

Ni–Ni coordination number CNNi–Ni as a function of time for (a) 6Hbulk/Ni(110) and (b) 6Hsurf/Ni(110). Each curve of the same temperature corresponds to one Ni atom, i.e., the number of markers at each time t is the total number of Ni nearest neighbors. (c) Distribution of Hsurf adsorption sites at 100 K, 200 K, and 400 K for 6Hsurf/Ni(110). (d) Examples of bridge (twofold coordinated) and hollow (threefold coordinated) sites for Hsurf shown with the orange and red bonds, respectively.

Close modal

To understand the effect of the unique properties of Hbulk/Ni(110) on adsorbates, we added a single CO2 molecule to the surface of the Hsurf/Ni(110) and Hbulk/Ni(110) systems. CO2 was placed in the well-known V-shaped chemisorbed structure on the short bridge-long bridge site [see Figs. 2(d)2(f)], which has been shown to be the most stable for CO2 chemisorption on Ni(110).21 We first examine how CO2 adsorption affects the Ni(110) surface structure, with rNi–Ni, CNNi–Ni, and Δd12 for each Hsurf/Ni(110) and Hbulk/Ni(110) system shown in Table II for the optimized geometries, as well as the two CO bond lengths rC–O,1 and rC–O,2. We see that CO2 adsorption has a negligible effect on rNi–Ni, with a maximum change of 0.07 Å, but we still observe an overall increase with H coverage relative to rNi–Ni of pristine Ni(110). In terms of Ni–Ni coordination, we see either no change or a reduction in CNNi–Ni toward the value for pristine Ni(110) (CNNi–Ni = 7). Interestingly, we see expansion of the interlayer spacing for all systems, which increases with coverage for both Hbulk and Hsurf. With respect to the CO bond lengths, we see that one of the CO bonds elongates, as typical of CO2 activated by metal surfaces,63 with increasing coverage for both Hbulk and Hsurf.

TABLE II.

Ni–Ni bond lengths rNi–Ni and Ni–Ni coordination numbers CNNi–Ni of surface Ni atoms and layer relaxation Δd12 with respect to the pure Ni(110) interlayer spacing.

rNi–Ni (Å)rC–O,1 (Å)rC–O,2 (Å)CNNi–NiΔd12 (%)
Ni(110) 2.45 1.25 1.27 7.0  
2Hsurf/Ni(110) 2.47 1.29 1.30 7.0 8.14 
2Hbulk/Ni(110) 2.48 1.30 1.31 7.0 7.95 
4Hsurf/Ni(110) 2.48 1.30 1.29 7.2 8.36 
4Hbulk/Ni(110) 2.51 1.39 1.32 7.0 13.7 
6Hsurf/Ni(110) 2.50 1.36 1.28 7.0 9.90 
6Hbulk/Ni(110) 2.51 1.35 1.28 7.0 10.75 
rNi–Ni (Å)rC–O,1 (Å)rC–O,2 (Å)CNNi–NiΔd12 (%)
Ni(110) 2.45 1.25 1.27 7.0  
2Hsurf/Ni(110) 2.47 1.29 1.30 7.0 8.14 
2Hbulk/Ni(110) 2.48 1.30 1.31 7.0 7.95 
4Hsurf/Ni(110) 2.48 1.30 1.29 7.2 8.36 
4Hbulk/Ni(110) 2.51 1.39 1.32 7.0 13.7 
6Hsurf/Ni(110) 2.50 1.36 1.28 7.0 9.90 
6Hbulk/Ni(110) 2.51 1.35 1.28 7.0 10.75 

Now, examining the electronic structure, in Table III, we show the charge on each H atom, qH, the Ni surface charge, qNi, and the charge of CO2, qCO2, as well as their changes relative to the H/Ni systems without CO2. It is well-known that chemisorbed CO2 on nickel is negatively charged, with qCO2 = −0.8–0.9e.21 According to Table III, our values for qCO2 are in excellent agreement with the literature. We observe that CO2 on Hbulk/Ni(110) is more negatively charged than CO2 on Hsurf/Ni(110) at all coverages, corresponding to increased charge transfer from the nickel surface to adsorbed CO2. This fits very well with our hypothesis from the previous section H/Ni that when H is in the bulk, the Ni(110) surface is more reductive than when H is on the surface. This is further substantiated by the computed work functions of 4Hbulk/Ni(110) and 6Hbulk/Ni(110), which are lower than those of 4Hsurf/Ni(110) and 6Hsurf/Ni(110). Due to the increased charge transfer from Ni to CO2, the charge of the first layer Ni atoms, qNi, of Hbulk/Ni(110) is greater than qNi of Hsurf/Ni(110) for all coverages, and Hbulk are slightly more negative than Hsurf. We also see that CO2 adsorption increases the charge on each Hsurf while decreasing the charge on each Hbulk by ∼0.02e per H atom.

TABLE III.

Charges per hydrogen atom qH, total Ni surface charge qNi, and charge on CO2, qCO2.

qHqNIqCO2
Ni(110)  0.25 −0.8 
2HSURF/NI(110) −0.21 (0.02) 1.20 (0.76) −0.83 
2HBULK/NI(110) −0.25 (−0.02) 1.01 (0.8) −0.82 
4HSURF/NI(110) −0.20 (0.02) 1.40 (0.64) −0.84 
4HBULK/NI(110) −0.23 (0.0) 1.12 (0.91) −0.97 
6HSURF/NI(110) −0.21 (0.02) 1.77 (0.63) −0.84 
6HBULK/NI(110) −0.24 (−0.02) 1.09 (0.84) −0.95 
qHqNIqCO2
Ni(110)  0.25 −0.8 
2HSURF/NI(110) −0.21 (0.02) 1.20 (0.76) −0.83 
2HBULK/NI(110) −0.25 (−0.02) 1.01 (0.8) −0.82 
4HSURF/NI(110) −0.20 (0.02) 1.40 (0.64) −0.84 
4HBULK/NI(110) −0.23 (0.0) 1.12 (0.91) −0.97 
6HSURF/NI(110) −0.21 (0.02) 1.77 (0.63) −0.84 
6HBULK/NI(110) −0.24 (−0.02) 1.09 (0.84) −0.95 

While the previous results give a static picture of CO2 reduction on H/Ni(110), experiments typically run at temperature well above 300 K.64 To examine temperature effects on the catalytic properties of these systems, we performed a series of AIMD simulations at 100, 400, and 1000 K for 6Hbulk/Ni(110) and 6Hsurf/Ni(110). As in the 400 K simulation for 6Hbulk/Ni(110) without CO2 (see Fig. S4a), we observe the emergence of single Hbulk to the surface on the same timescale, as shown in Fig. S4b. However, we also see that more Hbulk atoms have moved deeper into the bulk (layer 4), most likely due to the stabilization of CO2 on the surface and a lower number of available surface sites to which Hbulk can emerge, due to CO2. Thus, Hbulk and CO2 stabilize each other: Hbulk stabilizes CO2 by reducing the work function of the surface, enabling it to be more reduced than in the presence of Hsurf, and CO2 stabilizes Hbulk by taking up space on the surface. In fact, we observe diffusion of Hbulk deeper into the bulk (away from CO2) before any H emerges to the surface (see Fig. S4b). We also observe such diffusion deeper into the bulk even at low temperatures, as in the simulation at 100 K for 6Hbulk/Ni(110), but no emergence to the surface. We also estimated the lifetime of Hsub, which is very small, on the order of ∼1 ps, even in the 400 K simulation, confirming our hypothesis that the energy gained from moving from the bulk to the subsurface is enough to overcome the barrier to emerge to the surface.

At 1000 K, almost all Hbulk diffuse to the surface within 30 ps, as shown in Fig. 6(a). However, we observe that adsorption of CO2 also pushes one Hsurf atom to the bottom surface; see the comparison in Fig. 6. We attribute this to the reduction and stabilization of CO2, which dissociates within the first picosecond, as shown in Fig. 7. Interestingly, CO2 dissociated before any Hbulk diffused to the surface, i.e., in the absence of Hsurf and of Hsub, as shown in Fig. 7(b). The dissociation was subsequently followed by one Hbulk diffusing to the top surface and one Hbulk diffusing to the bottom surface, all during equilibration, and remained for the length of the trajectory. In contrast, when all the hydrogens start out as Hsurf, we observe that CO2 desorbs from the surface within the first picosecond of equilibration, as shown in Fig. S5. There is also a significant surface reconstruction for the CO2 + 6Hsurf/Ni(110) system at this temperature. This, coupled with the high work function of this system, destabilizes CO2 and causes it to desorb.

FIG. 6.

Height (z) of each H atom over time at 1000 K for 6Hbulk/Ni(110) with CO2 (a) and for 6Hbulk/Ni(110) without CO2 (b). Each curve in each plot corresponds to a single H atom. The black dashed lines correspond to the Ni layers, with the topmost black dashed line corresponding to the surface. Sample snapshots from the AIMD simulations at 1000 K are shown (c) with CO2 and (d) without CO2.

FIG. 6.

Height (z) of each H atom over time at 1000 K for 6Hbulk/Ni(110) with CO2 (a) and for 6Hbulk/Ni(110) without CO2 (b). Each curve in each plot corresponds to a single H atom. The black dashed lines correspond to the Ni layers, with the topmost black dashed line corresponding to the surface. Sample snapshots from the AIMD simulations at 1000 K are shown (c) with CO2 and (d) without CO2.

Close modal
FIG. 7.

(a) Side view of configuration at which CO2 dissociates during 1000 K simulation. (b) CO bond length, rCO (orange), that breaks during 1000 K equilibration run of 6Hbulk/Ni(110) and height z (green) of the first H atom to emerge to the surface, both as a function of time. The black dashed line corresponds to the top Ni surface.

FIG. 7.

(a) Side view of configuration at which CO2 dissociates during 1000 K simulation. (b) CO bond length, rCO (orange), that breaks during 1000 K equilibration run of 6Hbulk/Ni(110) and height z (green) of the first H atom to emerge to the surface, both as a function of time. The black dashed line corresponds to the top Ni surface.

Close modal

Schatz et al. showed that on bare Ni(110), CO2 dissociation and desorption are competing processes as the difference between desorption and dissociation energies is only 0.03 eV, while the presence of Hsub increases this difference to 0.13 eV, favoring dissociation.21 However, they claimed that Hbulk (any hydrogen below the second layer of Ni) does not alter the electronic structure of the surface in any significant way and thus does not affect adsorbed species. Our AIMD simulations offer an additional scenario that relates to the dynamic behavior of the system when temperature effects are taken into account. Our results show stabilization and reduction in CO2 in the presence of Hbulk, which has not been discussed in the literature previously. We also note that Schatz and others21,32,65 assumed that the benefits (e.g., higher reactivity) of Hsub are primarily due to the extra energy released when Hsub emerges to the surface prior to CO2 dissociation; however, our temperature-dependent simulations suggest that CO2 can dissociate before any emergence of Hsub/Hbulk as long as it can adsorb to the surface. To provide further evidence, we computed the adsorption energy of CO2 on optimized pure Ni(110) and 6H/Ni(110) systems, according to EadsCO2=ECO2+surfEsurf+ECO2. Here, ECO2+surf is the energy of the entire system consisting of CO2 adsorbed on pristine Ni(110) or one of the 6H/Ni systems, Esurf is the energy of the system without CO2, and ECO2 is the energy of a single CO2 molecule in the gas phase. For CO2 adsorption on pristine Ni(110), we obtain a value of −0.4 eV, which agrees well with the literature (−0.47 eV),21,32 indicating that at room temperature and above, there will be little to no adsorption of CO2 on the surface. For 6Hbulk/Ni(110), we obtain a stronger adsorption energy of −1.3 eV, while we obtain a positive value of 1.93 eV for 6Hsurf/Ni(110), suggesting that CO2 cannot adsorb on our 6Hsurf/Ni(110), leading to desorption from the surface in AIMD. Finally, we note that there is an additional physisorbed state with a very small adsorption energy (−0.1 eV) arising from dispersive van der Waals interactions, which is not likely to contribute to any relevant catalytic pathways at elevated temperatures.

In general, splitting CO2 into CO and O is thermodynamically unfavorable but becomes achievable in the presence of hydrogen, according to the following reverse water–gas shift reaction (rWGS):66 

CO2g+H2gCOg+H2Og,ΔH2980=9.8kcal/mol.

This reaction is known to occur through two mechanisms:21 (i) an associative mechanism, where CO2 is first hydrogenated to HCOO, which then dissociates into CO and OH, followed by further hydrogenation of OH to H2O, and (ii) a redox mechanism, where CO2 dissociates first into CO and O, followed by two subsequent hydrogenations of O to H2O. It has been demonstrated that the associative mechanism dominates in the absence of Hsub, while the redox mechanism, which has a lower energy barrier, dominates in the presence of Hsub.21 In our case, we see a preference for the redox mechanism in the presence of Hbulk as well. However, in previous studies, overcoming the CO2 dissociation barrier was attributed to the emergence to the surface of Hsub, but in our simulation, there is no emergence to the surface prior to CO2 dissociation. On pure Ni(110), temperature-programmed desorption (TPD) experiments show that chemisorbed CO2 is stable up to ∼220 K.67 At this temperature, CO2 partially desorbs and partially dissociates into CO and O. Because experiments typically run at temperatures well above 300 K,64 it is desirable to find ways to favor CO2 dissociation over CO2 desorption, particularly via the redox mechanism, as it has a lower energy barrier than the associative mechanism. We attribute the facile CO2 dissociation observed in our simulation to enhanced CO2 reduction and stabilization, which arises from the unique combination of several desirable properties, namely, a low work function and high d-band center. Therefore, we propose the possibility of a new redox mechanistic pathway for CO2 reduction, which may have a lower barrier to CO2 dissociation than previously identified pathways, as displayed in Fig. 8(a). A number of distinct 2Hbulk/Ni(110) (15) were optimized, and the most stable (No. 12) was used to generate six additional 4Hbulk/Ni(110) configurations that were also optimized; see new Fig. S9 and Table S3 of the supplementary material. Configuration #4 was used as the reference state for the reaction CO2(g) + H2(g) → H2O(g) and CO(g). We also calculated all the activation energy barriers that connect the intermediate products. As can be seen, the highest intermediate barriers are on the order of 1 eV, which are typically surmountable at temperature above 500 K. The largest energy penalty is the tandem desorption of the H2O and CO products, ∼ 3.3 eV. However, these events do not have to happen in tandem. Desorption of H2O alone amounts to about 0.8 eV, which is comparable to values in the literature (∼0.6 eV). For CO, we compute a desorption energy of about 2.7 eV, similar to other reported values (∼2.1 eV from the work of Schatz et al.). Both of these values are somewhat overestimated most likely due to the finite size basis employed or functional. In any event, CO desorption will require high temperatures (consistent with the conclusions of Schatz et al.), and as result, CO may remain on the surface for further hydrogenation. This reference allows us to also write the reaction pathway as a catalytic cycle; see Fig. 8(b). The main difference between our proposed pathway and previously identified pathways is the order of CO2 dissociation and H emergence to the surface.

FIG. 8.

Proposed redox mechanism for CO2 dissociation in the presence of Hbulk based on relative energies of the intermediates. (a) Potential energy landscape with barriers shown for selected steps. (b) Catalytic cycle representation of the reaction pathway.

FIG. 8.

Proposed redox mechanism for CO2 dissociation in the presence of Hbulk based on relative energies of the intermediates. (a) Potential energy landscape with barriers shown for selected steps. (b) Catalytic cycle representation of the reaction pathway.

Close modal

In the redox mechanism, the rate-determining step is CO2 dissociation, which is facile in our case due to the decreased work function of Hbulk/Ni(110). Following dissociation, two hydrogenations of O occur to give water; thus, it is required that two hydrogens emerge from the subsurface or bulk. We previously showed that the barrier for emergence from the bulk is comparable to that of emergence from the surface and on the order of 0.3 eV. As expected, in Fig. 8, the barrier for the first emergence is low, while the emergence of the second hydrogen is substantially increased, which was corroborated by Schatz et al.21 The first hydrogenation of O to hydroxyl and the second hydrogenation to water have been shown to be fast on Ni(110),68 with barriers of 0.57 and 1.03 eV, respectively,21 which correspond well with our computed barriers of 1.01 and 0.95 eV. As shown in Fig. 8, the only step (prior to desorption) that is significantly uphill is the final hydrogenation to water, Hsurf + CO∗ + OH∗ → CO∗ + H2O∗, as water is metastable on Ni(110). The last step of this pathway, desorption of the products CO and H2O, is also uphill, with a difference in the energy of 3.35 eV. Schatz et al. computed the desorption energies of CO and H2O on Hsub/Ni(110) to be 0.61 and 2.03 eV, respectively, both of which are larger than the corresponding desorption energies on pure Ni(110),21 which we attribute to the increase in adsorbate stability from Hsub.

This AIMD study demonstrates that bulk hydrogen plays a far more significant role in surface catalysis than previously depicted in the literature. In particular, we show that Hbulk can passively affect reaction mechanisms in new ways. This is beneficial for catalysis involving CO2 because keeping H inside (i) keeps the work function low, (ii) prevents accumulation of H on the surface, and (iii) enables the redox mechanism of rWGS, which has a lower reaction barrier than the associative mechanism. While the notion of bulk hydrogen influencing reactivity has been overlooked in the past because of its negligible effect on surface charge, we found that despite a similar surface charge to pristine Ni(110), the work function and d-band center of Hbulk/Ni(110) can be tuned to a desirable combination of low work function and high d-band center. While Hsurf/Ni(110) and Hsub/Ni(110) have higher d-band centers, their work functions also remain higher than those of Hbulk/Ni(110) at high coverages. We attribute this to differences in proximity to the Ni surface of each H species: moving the negatively charged H atoms deeper into the bulk, and thus away from the catalytic surface, reduces the work function with increasing H coverage, a trend which is not true for Hsurf and Hsub. This balance can greatly affect the behavior of adsorbed species, as we observed the dissociation of CO2 in the presence of Hbulk and the desorption of CO2 in the presence of Hsurf during molecular dynamics simulations. In addition, we also observed that CO2 dissociation can occur in the absence of Hsurf when only Hbulk is present, putting forth a reaction mechanism that has not been yet discussed in the catalysis community. Our results suggest that there is a much more complex interplay between Hsurf, Hsub, and Hbulk than previously anticipated. AIMD enabled the discovery of this newly found role of Hbulk and its potential to influence surface reactivity. This new reaction pathway of CO2 reduction could have implications for a variety of reactions involving CO2, such as methanol production, hydrogen production via the water–gas shift reaction, and CO2 consumption via the dry reforming reaction.

We dedicate this paper to the strong women in our lives, mothers, sisters, and daughters, that helped shape our work ethics and view challenges as opportunities. V.-A. G. had the incredible fortune to know the late Professor Pat Thiel during her Ph.D. studies at Iowa State and whose memory will be treasured.

See the supplementary material for the following: (1) Hartree potential of Ni(110); (2) comparison of Γ-point/k-point calculations; (3) charges per hydrogen atom qH, total Ni surface charge qNi, work function ϕ, and d-band center ϵd; (4) projected density of states of Ni d-orbitals; (5) MSD of H on the Ni(110)/6Hbulk system; (6) blue moon simulation details and energy decomposition; (7) H solubility in Ni; and (8) configurational ensemble for Ni(110)/Hbulk.

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, U.S. Department of Energy, Award/Contract Number 47319 (Multifunctional Catalysis to Synthesize and Utilize Energy C). The authors acknowledge computer resources through Research Computing (RI) located at Pacific Northwest National Laboratory (PNNL). PNNL is a multiprogram national laboratory operated for the DOE by Battelle under Contract No. DE-AC05-76RL01830. Computational resources were provided by a user proposal at the National Energy Research Scientific Computing Center located at Lawrence Berkeley National Laboratory. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

The data that support the findings of this study are available within the article and its supplementary material.

1.
R. J.
Kokes
and
P. H.
Emmett
,
J. Am. Chem. Soc.
81
,
5032
(
1959
).
2.
A. D.
Johnson
 et al,
Science
257
,
223
(
1992
).
3.
S. P.
Daley
 et al,
J. Am. Chem. Soc.
116
,
6001
(
1994
).
4.
K. L.
Haug
 et al,
J. Am. Chem. Soc.
120
,
8885
(
1998
).
5.
K.-A.
Son
and
J. L.
Gland
,
J. Am. Chem. Soc.
117
,
5415
(
1995
).
6.
A.
Stuck
 et al,
Phys. Rev. Lett.
74
,
578
(
1995
).
7.
M. E.
Dry
,
Appl. Catal., A
138
,
319
(
1996
).
8.
M. E.
Dry
,
J. Chem. Technol. Biotechnol.
77
,
43
(
2002
).
9.
W.
Wang
et al.,
Chem. Soc. Rev.
40
,
3703
(
2011
).
10.
S.
Kattel
 et al,
J. Am. Chem. Soc.
138
,
12440
(
2016
).
11.
S.
Kattel
,
P.
Liu
, and
J. G.
Chen
,
J. Am. Chem. Soc.
139
,
9739
(
2017
).
12.
P.
Munnik
,
P. E.
De Jongh
, and
K. P.
De Jong
,
J. Am. Chem. Soc.
136
,
7333
(
2014
).
13.
N. M.
West
et al.,
Coord. Chem. Rev.
255
,
881
(
2011
).
14.
S. T.
Ceyer
,
Acc. Chem. Res.
34
,
737
(
2001
).
15.
R. A.
Köppel
,
C.
Stöcker
, and
A.
Baiker
,
J. Catal.
179
,
515
(
1998
).
16.
A.
Baiker
,
Appl. Organomet. Chem.
14
,
751
(
2000
).
17.
J.
Wambach
,
A.
Baiker
, and
A.
Wokaun
,
Phys. Chem. Chem. Phys.
1
,
5071
(
1999
).
18.
O. V.
Krylov
and
A. Kh.
Mamedov
,
Russ. Chem. Rev.
64
,
877
(
1995
).
19.
J.
Li
,
E.
Croiset
, and
L.
Ricardez-Sandoval
,
J. Mol. Catal. A: Chem.
365
,
103
(
2012
).
20.
S. D.
Miller
and
J. R.
Kitchin
,
Surf. Sci.
603
,
794
(
2009
).
21.
W.
Lin
,
K. M.
Stocker
, and
G. C.
Schatz
,
J. Phys. Chem. C
120
,
23061
(
2016
).
22.
A.
Johnson
 et al,
Phys. Rev. Lett.
67
,
374
(
1991
).
23.
K. J.
Maynard
 et al,
Faraday Discuss.
91
,
437
(
1991
).
24.
C. N.
Satterfield
,
Heterogeneous Catalysis in Industrial Practice
, 2nd ed. (
McGraw Hill Book Co.
,
New York, NY
,
1991
).
25.
H. A.
Smith
,
A. J.
Chadwell
, and
S. S.
Kirslis
,
J. Phys. Chem.
59
,
820
(
1955
).
26.
J.-X.
Liu
 et al,
J. Am. Chem. Soc.
135
,
16284
(
2013
).
27.
W.
Chen
 et al,
ACS Catal.
7
,
8061
(
2017
).
28.
O. R.
Inderwildi
,
S. J.
Jenkins
, and
D. A.
King
,
J. Phys. Chem. C
112
,
1305
(
2008
).
29.
Y.
Men
 et al,
Catal. Today
125
,
81
(
2007
).
30.
E.
Vesselli
 et al,
J. Phys. Chem. C
115
,
1255
(
2011
).
31.
B.
Bhatia
and
D. S.
Sholl
,
J. Chem. Phys.
122
,
204707
(
2005
).
32.
G.
Peng
 et al,
Surf. Sci.
606
,
1050
(
2012
).
33.
A. P.
Ashwell
 et al,
J. Am. Chem. Soc.
139
,
17582
(
2017
).
34.
J. R.
Kitchin
 et al,
Phys. Rev. Lett.
93
,
156801
(
2004
).
35.
M.
Mavrikakis
,
B.
Hammer
, and
J. K.
Nørskov
,
Phys. Rev. Lett.
81
,
2819
(
1998
).
36.
F.
Solymosi
,
J. Mol. Catal.
65
,
337
(
1991
).
37.
S.-G.
Wang
 et al,
J. Phys. Chem. B
109
,
18956
(
2005
).
38.
A.
Mohsenzadeh
,
K.
Bolton
, and
T.
Richards
,
Surf. Sci.
627
,
1
(
2014
).
39.
A.
Michaelides
,
P.
Hu
, and
A.
Alavi
,
J. Chem. Phys.
111
,
1343
(
1999
).
40.
X.
Shen
 et al,
Phys. Chem. Chem. Phys.
19
,
3557
(
2017
).
41.
F.
Zhai
 et al,
J. Chem. Phys.
149
,
174704
(
2018
).
42.
G.
Peng
 et al,
J. Phys. Chem. C
116
,
3001
(
2012
).
43.
J.
Hutter
 et al,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
15
(
2014
).
44.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
,
Phys. Rev. Lett.
77
,
3865
(
1996
).
45.
S.
Grimme
 et al,
J. Chem. Phys.
132
,
154104
(
2010
).
46.
S.
Goedecker
,
M.
Teter
, and
J.
Hutter
,
Phys. Rev. B
54
,
1703
(
1996
).
47.
M.
Krack
,
Theor. Chem. Acc.
114
,
145
(
2005
).
48.
G.
Lippert
,
J.
Hutter
, and
M.
Parrinello
,
Mol. Phys.
92
,
477
(
1997
).
49.
J.
Vandevondele
 et al,
Comput. Phys. Commun.
167
,
103
(
2005
).
50.
D.
Frenkel
and
B.
Smit
,
Understanding Molecular Simulation: From Algorithms to Applications
(
Academic Press
,
2002
).
51.
S.
Nosé
,
J. Chem. Phys.
81
,
511
(
1984
).
52.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
53.
G.
Henkelman
,
A.
Arnaldsson
, and
H.
Jónsson
,
Comput. Mater. Sci.
36
,
354
(
2006
).
54.
K. E.
Laidig
and
R. F. W.
Bader
,
J. Chem. Phys.
93
,
7213
(
1990
).
55.
E.
Sanville
 et al,
J. Comput. Chem.
28
,
899
(
2007
).
56.
W.
Tang
,
E.
Sanville
, and
G.
Henkelman
,
J. Phys.: Condens. Matter
21
,
084204
(
2009
).
57.
M.
Yu
and
D. R.
Trinkle
,
J. Chem. Phys.
134
,
064111
(
2011
).
58.
B.
Cordero
 et al,
Dalton Trans.
2832
(
2008
).
59.
N.
Doudin
 et al,
ACS Catal.
9
,
7876
(
2019
).
60.
B.
Hammer
and
J. K.
Norskov
,
Nature
376
,
238
(
1995
).
61.
G.
Alefeld
and
J.
Völkl
,
Hydrogen in Metals I
(
Springer
,
Berlin
,
1978
), p.
321
.
62.
J.
Greeley
and
M.
Mavrikakis
,
Surf. Sci.
540
,
215
(
2003
).
63.
V.-A.
Glezakou
and
L. X.
Dang
,
B. P.
McGrail
,
J. Phys. Chem. C
113
,
3691
(
2009
).
64.
E.
Vesselli
 et al,
J. Am. Chem. Soc.
130
,
11417
(
2008
).
65.
V.
Ledentu
,
W.
Dong
, and
P.
Sautet
,
J. Am. Chem. Soc.
122
,
1796
(
2000
).
66.
A. D.
King
,
D. B.
Yang
, and
D.
Yang
,
J. Am. Chem. Soc.
102
,
1028
(
1980
).
67.
X.
Ding
 et al,
Phys. Rev. B
76
,
195425
(
2007
).
68.
E.
Vesselli
 et al,
J. Chem. Phys.
122
,
144710
(
2005
).

Supplementary Material