PNAS PLUS

Ab initio molecular dynamics of solvation effects on reactivity at electrified interfaces Jeffrey A. Herrona, Yoshitada Morikawab,1, and Manos Mavrikakisa,1 a Department of Chemical and Biological Engineering, University of Wisconsin–Madison, Madison, WI 53706; and bDivision of Precision Science & Technology and Applied Physics, Osaka University, Osaka 565-0871, Japan

Edited by Alexis T. Bell, University of California, Berkeley, CA, and approved June 29, 2016 (received for review March 20, 2016)

AIMD

| DFT | methanol | Pt(111) | electrocatalysis

M

ethanol presents a promising fuel alternative to hydrogen for low-temperature polymer electrolyte membrane fuel cells (1, 2). However, the slower reaction kinetics, compared with hydrogen, requires high overpotentials, which reduces the overall energy efficiency of these fuel cells. Additionally, CO molecules that are formed during the reaction poison Pt catalysts, limiting their activity (2). To address this problem, improved electrocatalysts must be developed (3, 4). Toward this goal, it is essential to understand the sequence of elementary steps that comprise the reaction mechanism. Through knowledge of the detailed mechanism on a particular catalyst, one can determine the ratedetermining step(s). In turn, this information can be used to design improved catalysts that are targeted toward facilitating the rate-determining step. First-principles density functional theory (DFT) calculations have become invaluable for working toward understanding these detailed reaction mechanisms at the atomic scale. In particular, these methods have successfully been used to calculate the thermochemistry (reaction energies) and kinetics (activation energy barriers) of a variety of different catalytic reactions in the vapor phase. However, modeling electrocatalytic reactions poses a greater challenge due to the complexity of the reaction environment at the active site. That is, the reaction occurs on a charged substrate that is surrounded by solvent (and electrolyte). A variety of methods have been developed to model the effects of solvent and charged surfaces on the thermochemistry and kinetics of electrocatalytic reactions (5–13), and the strengths of these approaches were summarized in a recent publication (13).

www.pnas.org/cgi/doi/10.1073/pnas.1604590113

Although these methods have made advances toward modeling this complicated active site, one of the many deficiencies is the inability to capture the effect of dynamic solvent reorganization (14). Here, we use ab initio (constrained) molecular dynamics (AIMD) to study the effect of water solvent reorganization and substrate charge toward methanol electrooxidation on a model Pt surface, Pt(111). In particular, we use a Blue Moon Ensemble (15– 20) to evaluate the free energy of activation for breaking the C−H and O−H bonds in methanol. This systematic study applies AIMD to a heterogeneously catalyzed reaction accounting for both solvent and charged surfaces. Results and Discussion In this section, we report the results of our AIMD simulation. The results are organized into three sections. First, we discuss calculated reaction thermochemistry and binding energies of adsorbates on uncharged, solvated Pt(111). Then, we present calculated activation energy barriers using the Blue Moon Ensemble for uncharged, solvated Pt(111). Finally, we present activation energy barriers for the same elementary steps on a charged, solvated Pt(111) surface using the Blue Moon Ensemble. Reaction Thermochemistry on Uncharged, Solvated Pt(111). First, we calculated the minimum energy binding geometries and the corresponding binding energies of methanol, methoxy (OCH3), and hydroxymethyl (CH2OH) on Pt(111) in the absence of additional charge or solvation. The equation for the binding energy (BE) was BE = Eadsorbate − Eclean − Eisolated, where Eadsorbate is the total energy of the adsorbate and surface when the adsorbate is adsorbed, Eclean is the total energy of the clean slab, and Eisolated is the total energy of the isolated gas-phase adsorbate. We found

Significance Low-temperature fuel cells are efficient energy conversion devices that face a number of hurdles toward commercialization, including difficulties in storing hydrogen. Methanol represents a liquid-phase fuel alternative to hydrogen, yet the high cost of Pt-based catalysts limits fuel cells’ economic viability. Toward improved, lower-cost catalyst design, a fundamental understanding of the methanol electrooxidation reaction mechanism is necessary. Density functional theory calculations have become invaluable in elucidating these reaction mechanisms, although the complex reaction environment including solvation of a charged electrode has been a challenge to model. Using ab initio molecular dynamics, via the Blue Moon Ensemble, we have investigated methanol electrooxidation on a solvated and charged Pt(111) surface to understand the effect of solvation and charge on the reaction energetics. Author contributions: Y.M. and M.M. designed research; J.A.H. performed research; J.A.H., Y.M., and M.M. analyzed data; and J.A.H., Y.M., and M.M. wrote the paper. The authors declare no conflict of interest. This article is a PNAS Direct Submission. 1

To whom correspondence may be addressed. Email: [email protected] or morikawa@ prec.eng.osaka-u.ac.jp.

PNAS | Published online August 8, 2016 | E4937–E4945

ENGINEERING

Using ab initio molecular dynamics as implemented in periodic, self-consistent (generalized gradient approximation Perdew– Burke–Ernzerhof) density functional theory, we investigated the mechanism of methanol electrooxidation on Pt(111). We investigated the role of water solvation and electrode potential on the energetics of the first proton transfer step, methanol electrooxidation to methoxy (CH3O) or hydroxymethyl (CH2OH). The results show that solvation weakens the adsorption of methoxy to uncharged Pt(111), whereas the binding energies of methanol and hydroxymethyl are not significantly affected. The free energies of activation for breaking the C−H and O−H bonds in methanol were calculated through a Blue Moon Ensemble using constrained ab initio molecular dynamics. Calculated barriers for these elementary steps on unsolvated, uncharged Pt(111) are similar to results for climbing-image nudged elastic band calculations from the literature. Water solvation reduces the barriers for both C−H and O−H bond activation steps with respect to their vapor-phase values, although the effect is more pronounced for C−H bond activation, due to less disruption of the hydrogen bond network. The calculated activation energy barriers show that breaking the C−H bond of methanol is more facile than the O−H bond on solvated negatively biased or uncharged Pt(111). However, with positive bias, O−H bond activation is enhanced, becoming slightly more facile than C−H bond activation.

Table 1. Average (GGA-PBE) binding energies () of CH3OH, CH2OH, and OCH3 on solvated, uncharged Pt(111) derived from AIMD simulation Vapor phase

Parameter CH3OH, eV CH2OH, eV OCH3, eV CH3OH → CH2OH + H CH3OH → OCH3 + H

This work

Literature

Aqueous phase, uncharged

−0.09 −2.10 −1.58 −0.04 0.76

−0.33*, −0.61† −1.98*, −1.93† −1.54* −0.16* 0.62*, 0.59†

−0.14 −2.09 −1.00 0.02 1.53

Included is thermochemistry () for CH3OH dehydrogenation to CH2OH and OCH3 on uncharged Pt(111) with hydrogen balance from isolated H2(g). Vapor-phase binding energies on uncharged Pt(111) are presented along with previous theoretical and experimental results for comparison. *Perdew–Wang 91 (PW91) GGA DFT calculations from refs. 23 and 24. † Microcalorimetry experiments from ref. 22.

generalized gradient approximation Perdew–Burke–Ernzerhof (21) (GGA-PBE) binding energies of −0.09 eV (methanol, physisorbed), −1.58 eV (OCH3, top site), and −2.10 eV (CH2OH, top site), for the minimum energy structure of each species. Note that we found that CH2OH is more stable than OCH3 in the gas phase by 0.27 eV, and it is more stable by 0.80 eV on Pt(111), in the absence of solvation. These values are compared with experimental (22) and theoretical (23, 24) results from the literature in Table 1. We then explored how solvation in water affects these binding energies. To do this, we performed unconstrained AIMD simulations of the adsorbate over the Pt(111) surface that is surrounded by four layers of water. To maintain the density in the aqueous phase close to 1 g/cm3, we removed a single water molecule (corresponding to 15 water molecules) when adding an adsorbate. We calculated the average energy of the slab with an adsorbate and 15 water molecules (hEadsorbate i), and the slab with 16 water molecules (hEclean i). The average binding energy (hB. E. i) was calculated

from the following equation: hB. E. i = hEadsorbate i − hEclean i − Eisolated + μH2 O, where μH2 O is the chemical potential of water, and Eisolated is the total energy of the isolated gas-phase adsorbate. The chemical potential of water was calculated from bulk water via AIMD (2-ps trajectory) using a 9.85 Å × 9.85 Å × 9.85 Å cubic supercell with 32 water molecules (1.00 g/cm3 density). For reference, from our simulation we derived a chemical potential for bulk water of −0.67 eV. The ultimate goal of our study is to understand how solvation and surface charge affect adsorption and reactions on the surface. Therefore, we tested the convergence of the average energy, with respect to averaging time length, using CH2OH* with 15 water molecules on Pt(111) with 0.15 e− in excess charge in the surface (Fig. 1). Previous studies have found that, by adding or subtracting, this electron density corresponds to a potential range of ±0.4 VSHE (25). This calculation was chosen, rather than a clean surface with 16 H2O molecules, to ensure that both the adsorbate reorganization on the surface and water dynamics were converged. From the results, we see that, as expected, the oscillations in the average energy are dampened as the averaging time length is increased. Averaging the total energy over a 2-ps simulation (after allowing the simulation to equilibrate, i.e., locate the local minima) was deemed sufficient because the calculated average energy fluctuated by plus or minus ∼0.1 eV around the more stable average obtained with 3-ps averaging time. For other calculations, we used the equilibrated water geometry obtained here as an initial guess to reduce the equilibration time. Following this procedure, we calculated the average binding energies of CH3OH, CH2OH, and OCH3 as shown in Table 1. The corresponding AIMD trajectories are shown in Fig. 2. From these results, we see that the binding energies of CH3OH and CH2OH on uncharged Pt(111), in the aqueous phase, are similar to the respective ones in the vapor-phase case. On the other hand, OCH3 is significantly destabilized (by 0.58 eV) in the aqueous phase. These differences are mainly due to the binding geometries and chemical functionalities of these adsorbates. In Fig. 3A, we show the optimized geometries of CH3OH, OCH3, and CH2OH on uncharged Pt(111), and in

Fig. 1. Convergence of average energy of CH2OH* (hECH2 OH i) with 15 H2O molecules on Pt(111) with 0.15 e− excess charge (per six surface Pt atoms) included. Running average of total energy (hECH2 OH i) over a varying length of averaging time, from 1 ps to 3.5 ps (see figure key) is shown. The time coordinate for the average energy corresponds to the point in simulation time from which the averaging window begins (e.g., for an averaging length of 1 ps, the average that is shown at 500 fs corresponds to the average from 500 fs to 1,500 fs). The total energy from AIMD simulation is shown in the dotted line.

E4938 | www.pnas.org/cgi/doi/10.1073/pnas.1604590113

Herron et al.

Fig. 3B, we show representative snapshots from the AIMD trajectories of the respective species. In the vapor phase, CH3OH binds with the CH3 group oriented toward the Pt(111) surface with the C–O bond axis 56° from the surface normal. The O–H bond is pointed toward the surface. Interestingly, in the aqueous phase, the C–O bond axis is nearly parallel with the surface, and the O–H bond is oriented toward the surface, rather than toward

Fig. 3. (A) Vapor phase—the optimized binding geometries of CH3OH, CH2OH, and OCH3 on uncharged Pt(111). (B) Aqueous phase—representative snapshots from AIMD trajectory of equilibrated CH3OH, CH2OH, and OCH3 on uncharged Pt(111). The snapshots were constructed by repeating the unit cell twice each in the x and y directions (2 × 2) to enhance clarity of the water network. The adsorbates are highlighted in red ellipses in the aqueous phase.

Herron et al.

PNAS | Published online August 8, 2016 | E4939

PNAS PLUS ENGINEERING

Fig. 2. Dynamic change of binding energies for CH3OH (black line), CH2OH (red line), and OCH3 (blue line) on uncharged solvated Pt(111). Average binding energies (hB. E. i) (averaged over 2 ps) of these adsorbates are shown with dotted lines.

the aqueous phase. A neighboring water molecule has its O–H bond pointing toward the O atom of CH3OH, forming a hydrogen bond. In the vapor phase, CH2OH is bound to a top site of Pt(111) through its C atom, the C–O bond axis is 68° away from the surface normal, and the O–H bond is oriented parallel to the surface. In the aqueous phase, the geometry is not significantly changed. The O–H bond is oriented slightly away from the surface, toward a neighboring water molecule, and the O–H bond of a neighboring water molecule is pointed toward the O atom, like with CH3OH. In the vapor phase, OCH3 is bound to the top site on Pt(111) through its O atom and the C–O bond oriented 57° with respect to the surface normal, exposing the CH3 group to the environment. In the aqueous phase, the equilibrated binding geometry is on the top site, where the C–O bond is oriented nearly parallel to the surface, like in the vapor phase. However, unlike CH2OH, there are no hydrogen-bonded water species. Therefore, the differences in binding energy trends due to water solvation are because of stabilizing hydrogen bonding in CH3OH and CH2OH, whereas OCH3 is not stabilized. From the calculated binding energies, we evaluated the thermochemistry,hΔEi, for dehydrogenation of CH3OH to CH2OH and to OCH3. The hydrogen balance was taken from vaporphase H2, i.e., CH3OH → CH2OH + 1/2H2 and CH3OH → OCH3 + 1/2H2. Therefore, the thermochemistry includes only solvation for methanol, methoxy, and hydroxymethyl species. For example, the hΔEi for methanol dehydrogenation to hydoxymethyl is given by hΔEi = hB. E. CH2 OH i − hB. E. CH3 OH i + Eisolated  CH2 OH + ð1=2ÞEisolated  H2 − Eisolated  CH3 OH , where Eisolated CH2OH, Eisolated CH3OH, and Eisolated H2 are the energies of isolated vapor-phase CH2OH, CH3OH, and H2, respectively. Substituting our definition of binding energy into this equation, the overall result is hΔEi = hECH2 OH i − hECH3 OH i + ð1=2ÞEisolated  H2 . As a result, the

Fig. 4. Constrained AIMD trajectories for CH3OH dehydrogenation to CH2OH on Pt(111) at (A) −0.4 VSHE with solvation, (B) 0.0 VSHE with solvation, (C) 0.4 VSHE with solvation, and (D) uncharged and without solvation. (E) The mean constraint force derived from the trajectories is plotted versus the reaction coordinate, C–H bond stretching.

only energies included are the average total energies taken from the AIMD simulations and the energy of isolated H2. The results are shown in Table 1. The reaction energies for CH3OH dehydrogenation to CH2OH on solvated, uncharged Pt(111) are similar to the vapor-phase energetics because the binding energies of CH3OH and CH2OH are similar in the two phases. However, dehydrogenation to OCH3 is significantly more endothermic in the aqueous phase, by 0.77 eV, due to the lack of stabilization of surface-adsorbed OCH3 by the aqueous phase. Free Energies of Activation for Vapor-Phase Chemistry on Clean Pt(111). First, we used constrained molecular dynamics to

evaluate the free energies of activation for CH 3 OH dehydrogenation on uncharged, Pt(111) without solvation using a Blue Moon Ensemble (15–20) as described in detail in Methods. This is analogous to vapor-phase methanol dehydrogenation. The chosen reaction coordinate involves stretching of the C–H (to form CH2OH) or O–H (to form OCH3) bonds. The purpose of this calculation is to compare the E4940 | www.pnas.org/cgi/doi/10.1073/pnas.1604590113

results versus the traditional nudged elastic band (NEB) method (26) and to isolate the effects of solvation, charge, and calculation method. Table 2. Free energies of activation () for CH3OH dehydrogenation to CH2OH and OCH3 on uncharged, solvated Pt(111) derived from Blue Moon Ensemble via constraining C–H and O–H bond lengths, respectively Vapor phase

Parameter CH3OH → CH2OH + H CH3OH → OCH3 + H

This work

Literature

Aqueous phase, uncharged

0.60 0.71

0.67* 0.81*

0.52 0.67

Vapor-phase results on unsolvated Pt(111) are calculated using the Blue Moon Ensemble for comparison purposes. *PW91 GGA DFT calculations from refs. 23 and 24. Note that product is adsorbed H rather than gas-phase H2.

Herron et al.

PNAS PLUS ENGINEERING

Fig. 5. Constrained AIMD trajectories for CH3OH dehydrogenation to OCH3 on Pt(111) at (A) −0.4 VSHE with solvation, (B) 0.0 VSHE with solvation, (C) 0.4 VSHE with solvation, and (D) uncharged and without solvation. (E) The mean constraint force derived from the trajectories is plotted versus the reaction coordinate, O–H bond stretching.

The constraint forces obtained from these AIMD simulations as a function of bond length are shown in Figs. 4 and 5 for C–H bond and O–H bond activation, respectively [we note that these figures also show the corresponding trajectories for the simulations using solvated and charged Pt(111), as will be described in Free Energies of Activation on Uncharged, Solvated Pt(111) and Free Energies of Activation on Charged, Solvated Pt(111). The corresponding activation energies are included in Table 2. CH3OH → CH2OH + H: Uncharged Pt(111) without solvation. The reaction coordinate resembles that of previous vapor-phase results using more traditional NEB methods (26), although the transition state occurs at a C–H bond length of 1.48 Å, whereas a transition state bond length of 1.56 Å has been calculated in the vapor phase using the NEB method (23, 24). We find a free energy of activation of 0.60 eV (at 300 K), whereas the NEB method predicts an activation energy barrier of 0.67 eV (at 0 K) (23, 24). This shows that one is able to accurately calculate activation energy barriers for surface phenomena using the Blue Moon Ensemble. CH3OH → OCH3 + H: Uncharged Pt(111) without solvation. As with C–H bond activation, the O–H bond activation reaction coordinate reHerron et al.

sembles that of previous vapor-phase results using NEB method (23, 24). However, the transition state O–H bond length is 0.18 Å shorter than what has been reported using the NEB method (23, 24). The vapor-phase barrier from the NEB method has been calculated as 0.81 eV (23, 24), 0.10 eV more than our calculated vapor-phase barrier using the AIMD approach. Free Energies of Activation on Uncharged, Solvated Pt(111). We used constrained molecular dynamics to evaluate the free energies of activation for CH3OH dehydrogenation on charged, solvated Pt(111) using a Blue Moon Ensemble (15–20), as described in detail in Methods. The chosen reaction coordinate involves stretching of the C–H (to form CH2OH) or O–H (to form OCH3) bonds; characteristic snapshots of the molecular geometry from each Blue Moon simulation are shown in Figs. 6 and 7, respectively. CH3OH → CH2OH + H: Solvated, uncharged Pt(111). When the C–H bond has been stretched from its equilibrium bond length of 1.06 Å to 1.22 Å, the carbon atom interacts weakly with the Pt(111) surface. (Characteristic snapshots from the constrained AIMD trajectories are shown in Fig. 6.) The CH3OH molecule remains close to the surface, while the orientation of the C–H PNAS | Published online August 8, 2016 | E4941

Fig. 6. Characteristic snapshots from constrained AIMD simulations of H–CH2OH bond breaking reaction on uncharged, solvated Pt(111) at fixed C–H bond lengths of 1.22 Å, 1.37 Å, and 1.53 Å. The snapshots were constructed by repeating the unit cell twice each in the x and y directions (2 × 2) to enhance clarity of the water network. Each snapshot shows approximately one unit cell length in the x direction (three Pt atoms wide). A second methanol molecule in each image was removed to improve clarity. The orientation of the camera with respect to the surface is approximately the same for each of these images. The stretched H–CH2OH molecule is highlighted in a red circle in each image, with the constrained H highlighted in yellow.

bond fluctuates considerably, that is, there is near-free rotation of the C–H bond axis. As the bond elongates to 1.37 Å, the C atom is positioned above a top site of Pt(111) with the elongated C–H bond oriented toward the top site. The C–O bond axis fluctuates from nearly parallel to the surface to almost perpendicular with the O–H bond hydrogen-bonded with the neighboring water molecules. When the C–H bond has been stretched to 1.53 Å, we are beyond the transition state. The C atom of CH2OH is bound to the top site of Pt(111) with the elongated C–H bond pointed toward a neighboring bridge site, where the H atom is bound. Based on linear interpolation of the calculated mean constraint force, the transition state (where the mean constraint force crosses from negative to positive; Fig. 4E) corresponds to a C–H bond length of 1.46 Å and an activation energy barrier of 0.52 eV. The reaction coordinate resembles that of previous vaporphase results using more traditional NEB methods (26), although the transition state in the aqueous phase occurs at a C–H bond length that is ∼0.10 Å shorter (at 1.56 Å in the vapor phase using the NEB method) and with a barrier that is 0.15 eV less activated (24). Our Blue Moon Ensemble calculations on a Pt(111) surface without any solvation identified a transition state bond length of 1.48 Å. Therefore, the differences in bond lengths between the aqueous phase Blue Moon results and the vapor-phase NEB

results are likely due to differences in the specific computational methodologies rather than being an effect of solvation. Comparing the Blue Moon Ensemble vapor-phase result with the solvated result, we conclude that the activation energy barrier for H–CH2OH activation is reduced by 0.08 eV as a result of aqueous solvation. CH3OH → OCH3 + H: Solvated, uncharged Pt(111). The equilibrium O–H bond length in CH3OH is 0.95 Å. As the bond is stretched slightly, to 1.11 Å, the molecule remains close to the Pt(111) surface, the O–H bond is oriented away from the surface, and hydrogen is bonded with water. (Characteristic snapshots from the constrained AIMD trajectories are shown in Fig. 7.) This continues up until the O–H bond is stretched to 1.43 Å, where the O atom becomes bound to the top site of Pt(111). Here, the O–H bond is oriented toward a neighboring top site. As the bond elongates further, the O atom binds to the top site with the O–H bond oriented toward a neighboring top site, although the H atom becomes bound to that neighboring top site rather than to the bridge site between the two top sites. Through linear interpolation of the mean constraint force, we estimate that the transition state corresponds to an O–H bond length of 1.58 Å with an activation energy barrier of 0.67 eV. Again, the reaction coordinate resembles that of previous vaporphase results using NEB methods (26), although the transition

Fig. 7. Characteristic snapshots from constrained AIMD simulations of CH3O–H bond-breaking reaction on uncharged, solvated Pt(111) at fixed O–H bond lengths of 1.11 Å, 1.27 Å, 1.43 Å, and 1.59 Å. The snapshots were constructed by repeating the unit cell twice each in the x and y directions (2 × 2) to enhance clarity of the water network. Each snapshot shows approximately one unit cell length in the x direction (3 Pt atoms wide). The second methanol molecule in each image was removed to improve clarity. The angle of orientation of the camera with respect to the surface normal in the latter three images is slightly increased with respect to the first image. The stretched CH 3O–H molecule is highlighted in a red circle in each image, with the constrained H highlighted in yellow.

E4942 | www.pnas.org/cgi/doi/10.1073/pnas.1604590113

Herron et al.

Solvation Center of Mass. To provide additional quantification to this solvation effect, we calculated the time-averaged, z-coordinate center of mass of oxygen atoms from the surrounding water molecules from the constrained AIMD simulations. The results are referenced relative to the z coordinate of the top layer of Pt atoms, as shown in Fig. 8. Additionally, we calculated the result from the surface without a methanol molecule in the unit cell and found a value of 8.40 Å. We note that, to maintain near-constant density in our simulations, when adding the methanol molecule, we removed one water molecule from our unit cell (leaving 15 water molecules, corresponding to approximately four layers). Therefore, the exact numerical values in this analysis should be read with caution. The z coordinate of the oxygen atoms’ center of mass provides an indication of how far the water molecules are from the surface or surface-bound adsorbate. For adsorbates that present a hydrophobic group to the surrounding water, we would expect the z coordinate to increase (i.e., the water molecules move away from the surface). For adsorbates that present a hydrophilic group to the surrounding water, we would expect a decrease in the z coordinate as the water molecules solvate the adsorbate. For O–H bond activation, as the O–H bond is stretched from 1.11 Å to 1.27 Å, the z coordinate of the oxygen center of mass increases from 8.34 Å to 9.17 Å as the water molecules move away from the surface-bound, hydrophobic methyl group. The z coordinate of the center of mass increases slightly more, by 0.12 Å, as the O–H bond is stretched past the transition state. For comparison, we calculated an equilibrium value of 9.89 Å from the unconstrained AIMD simulations of methoxy adsorbed on solvated uncharged Pt(111) (note that there is one fewer hydrogen atom in the binding energy simulations). For C–H bond activation, as the C–H bond is stretched to 1.37 Å, the z coordinate increases to 9.00 Å. Then, as the bond is stretched beyond the transition state (1.46 Å), the CH2 group is bound closer to the surface and the center of mass of the oxygen atoms moves closer to the surface, at 8.62 Å. From our unconstrained AIMD simulations of CH2OH adsorbed on solvated uncharged Pt(111), an equilibrium center of mass coordinate of 8.65 Å is found. We note that this value is much closer to our constrained AIMD result than the CH3O–H activation result (although we have not extended our constrained AIMD simulations past the transition state all of the way to the local minima final states, as they are in the unconstrained simulations). Importantly, these results show that the center of mass of oxygen is closer to the surface-bound CH2OH species than to the OCH3 species, except at very small stretching of the constrained bond. Based on linear interpolation of the results presented in Fig. 8, at the transition states (1.46 Å for H–CH2OH bond activation and 1.58 Å for CH3O–H bond activation), the center of mass of oxygen is much closer to the surface for C–H bond activation than O–H bond activation. Herron et al.

PNAS PLUS

Fig. 8. Time-averaged z-coordinate center of mass of oxygen atoms from water molecules (15 water molecules per unit cell) with respect to the Pt surface as a function of constrained bond length, calculated from constrained AIMD simulations on uncharged Pt(111). Black (squares) shows the results from stretching the O–H bond, and red (circles) shows the results from stretching the C–H bond. As the z coordinate of the center of mass increases, water molecules move away from the surface and the surface-bound adsorbate, corresponding to a weakening of hydrogen bond stabilization of the adsorbate. For reference, the transition state bond length is 1.46 Å for H–CH2OH and 1.58 Å for CH3O–H. The green dashed line shows the oxygen center of mass without methanol in the simulation (corresponding to 16 water molecules in the unit cell).

Overall, we find that breaking the C–H bond is more facile than breaking the O–H bond in methanol on solvated uncharged Pt(111). The transition state for breaking the C–H bond occurs after the C–H bond is stretched by 0.39 Å, whereas the O–H bond must be stretched by 0.50 Å. The surrounding aqueous environment stabilizes the reaction coordinates for both of these elementary steps. The barrier for C–H bond activation is 0.08 eV less with solvation, whereas O–H bond activation is only 0.04 eV less with solvation; the reason for this is that the OH group of CH2OH remains hydrogen-bonded throughout the elementary step. In contrast, the water stabilization of OCH3 diminishes once the O–H bond is stretched and O binds to the Pt surface. Our calculated oxygen center of mass for the surrounding water molecules relative to the Pt demonstrates that the water molecules move farther away from the surface as the O–H bond is increasingly stretched, whereas the water molecules do move away during C–H stretching, but only at intermediate stretching levels, coming back closer to the surface after the C–H bond in CH3OH is broken. Free Energies of Activation on Charged, Solvated Pt(111). Next, we used constrained molecular dynamics to evaluate the free energies of activation for CH3OH dehydrogenation on charged, solvated Pt(111) using a Blue Moon Ensemble. The Pt(111) surface was charged by adding or subtracting 0.15 e− to the surface. Previous studies have found that, by adding (or subtracting), this electron density corresponds to a potential range of −0.4 VSHE (+0.4 VSHE when removing electron density) (25). The chosen reaction coordinate involves stretching of the C–H (to form CH2OH) or O–H (to form OCH3) bonds. The reaction coordinates for the two elementary steps studied, with a charged surface, did not change significantly from the uncharged, solvated results (Figs. 6 and 7). The constraint forces obtained from these constrained AIMD simulations PNAS | Published online August 8, 2016 | E4943

ENGINEERING

state occurs at a bond length that is ∼0.2 Å shorter (1.80 Å in the vapor phase from the NEB method). Our vapor-phase result using the Blue Moon Ensemble predicts a transition-state bond length that is only 0.02 Å longer than the aqueous phase result. Therefore, comparing the vapor-phase result calculated using the Blue Moon Ensemble versus the solvated result, we find that solvation has little effect on the transition state bond length. Additionally, we find that solvation lowers the free energy of activation by only 0.04 eV (while the barrier for C–H bond activation was lowered by 0.08 eV). The preferential lowering of the C–H bond activation step suggests that water solvation tends to enhance selectivity toward the CH 2 OH deprotonation product (as opposed to the CH3O deprotonation product). As the O–H bond is elongated, it becomes oriented toward the surface, rather than the solvent, breaking the hydrogenbond interactions.

as a function of bond length (reaction coordinate) are shown in Figs. 4 and 5, respectively, and the corresponding activation energies are provided in Table 3. We note that the value at 0.0 VSHE corresponds to the uncharged results that were presented in Free Energies of Activation on Uncharged, Solvated Pt(111). CH3OH → CH2OH + H: Charged, solvated Pt(111). We find that the free energy of activating the C–H bond in CH3OH to form CH2OH is not significantly affected by the electrode potential. The barrier is ∼0.50 eV for all potentials probed, and the reaction coordinate is similar for the entire potential range we studied. Based on linear interpolation of the calculated mean constraint force, the transition state (where the mean constraint force crosses from negative to positive; Fig. 4E) corresponds to a C–H bond length of 1.49 Å at −0.4 V SHE, 1.46 Å at 0.0 VSHE, and 1.46 Å at +0.4 V SHE. We note that these differences are likely within the error of our calculation. CH3OH→OCH3 + H: Charged, solvated Pt(111). Unlike the H–CH2OH bond activation, here we find that the activation energy for CH3O–H bond activation is lowered at positive potentials. The calculated activation energy barrier is 0.62 eV at −0.4 VSHE and 0.67 eV at 0.0 VSHE, and decreases to 0.48 eV at +0.4 VSHE. Through linear interpolation of the mean constraint force, we estimate that the transition state corresponds to an O–H bond length of 1.58 Å at −0.4 VSHE, 1.58 Å at 0.0 VSHE, and 1.59 Å at +0.4 VSHE. We note that these bond length differences are likely within the error of our calculation. Therefore, although the reaction coordinate does not change noticeably (nor the transition state bond length), there is a ∼0.2-eV reduction in the activation energy barrier at positive potentials. Interestingly, we find that on uncharged and negatively biased Pt(111), C–H bond activation is more facile than the O–H bond activation. However, at +0.4 V, the barrier for O–H bond activation becomes slightly lower than that for C–H bond activation. We note that when either the C–H or O–H bonds are oriented away from the surface and are constrained at a significant stretch from their equilibrium, often a neighboring water molecule will spontaneously dissociate during the AIMD simulation. The H atom of the water binds to the strained methanol species, regenerating the C–H or O–H bond. At the same time, the constrained H atom binds to the dissociated water, regenerating its O–H bond. The end result is the formation of an unphysical, constrained methanol–water couple. This phenomena was found for both uncharged and charged Pt(111). As a result, we were unable to obtain an activation energy barrier for dissociation of these C–H or O–H bonds into the liquid phase directly (forming an H3O species). Rather, we were only able to obtain results where the C–H or O–H bonds break to form a surface-bound H; this may be the reason why we did not find a strong potential dependence for these barriers. Table 3. Free energies of activation () for CH3OH dehydrogenation to CH2OH and OCH3 on charged, solvated Pt(111) derived from Blue Moon Ensemble via constraining C–H and O–H bond lengths, respectively Aqueous phase Parameter CH3OH → CH2OH + H CH3OH → OCH3 + H

−0.4 V

0.0 V

0.4 V

0.48 0.62

0.52 0.67

0.51 0.48

Pt(111) surface was charged by adding or subtracting 0.15 e− to the surface. Previous studies have found that, by adding (or subtracting), this electron density corresponds to a potential range of −0.4 VSHE (+0.4 VSHE when removing electron density) (34).

E4944 | www.pnas.org/cgi/doi/10.1073/pnas.1604590113

Conclusions Using AIMD simulations, we investigated the first elementary steps for methanol electrooxidation in an aqueous environment on Pt(111). In particular, we compared the effect of electrode potential and solvation on the first proton/electron transfer event, electrooxidation of methanol to methoxy or hydroxymethyl. We found that for an uncharged, solvated Pt(111) surface, the binding energies of methanol and hydroxymethyl are similar to their vapor-phase values, but the binding of methoxy is destabilized by ∼0.6 eV. As a result, the reaction thermochemistry for methanol electrooxidation to methoxy is more endothermic in the aqueous phase, whereas electrooxidation to hydroxylmethyl is approximately isoenergetic when comparing the aqueous phase and vapor-phase results. Further, we calculated the activation energy barriers for breaking the C–H and O–H bond of methanol using the constrained molecular dynamics approach through a Blue Moon Ensemble in the vapor phase and in the aqueous phase [on both charged and uncharged Pt(111)]. The vapor phase results calculated using the Blue Moon Ensemble provide activation energy barriers that are similar to traditional static NEB methods. Furthermore, we find that solvation on uncharged Pt(111) lowers the barrier for C–H bond activation by 0.08 eV, while only lowering the barrier for O–H bond activation by 0.04 eV. For both the vapor-phase and aqueous phase results, C–H bond activation is more facile. When charging the surface, we find that C–H bond activation remains more facile than O–H bond activation with negative bias, but, with positive bias, O–H bond activation becomes slightly more facile. Methods All calculations in this work were performed using planewave DFT as implemented in the Simulation Tool for Atom Technology (STATE) (27, 28). The Pt(111) surface is represented using a periodic 3 × 2 unit cell with three atomic layers. The optimized Pt lattice constant of 3.949 Å was used, which is in close agreement with the experimental value (29), 3.92 Å. The ionic cores were described using ultrasoft Vanderbilt pseudopotentials (30, 31). The wave functions and the augmentation charge were expanded by a plane wave basis set with the cutoff energies of 25 Ry and 225 Ry, respectively. GGA-PBE was used for the exchange–correlation energy functional. The surface Brillouin zone was sampled with a 2 × 4 × 1 Monkhorst–Pack k-point mesh (32). The first-order Methfessel–Paxton scheme was used to deal with the Fermi level (33). The smearing width was set to 0.054 eV. To simulate the charged surface with solvation, the “effective screening medium” method was applied (14). In the surface-normal direction (in the z axis), a semiinfinite continuum with an infinite dielectric constant, i.e., a classical conductor, was located at z1 = ∼18.4 Å away from the center of the Pt atoms in the top layer of the surface (z = 0 Å), while another region (z < z1) was characterized by the dielectric constant of unity, i.e., the vacuum medium. In this region, we introduced 16 water molecules. To apply a potential bias, excess electrons per supercell were introduced, whereby the same amount of charge with the opposite sign was induced at the surface of the classical conductor. AIMD simulations were performed using STATE software (27, 28). The top two layers of surface atoms, along with all water molecules and adsorbates, were freely allowed to move, whereas the bottom layer of the slab was fixed. An artificial boundary was placed above the surface and water molecules, restricting their movement and maintaining the density at ∼1 g/cm3. In the molecular dynamics simulations, we have adopted the deuterium mass for H atoms to implement a longer time step (1.2 fs). Velocities were scaled to maintain the temperature around 300 K. Constrained molecular dynamics were performed while constraining the C–H or O–H bond in methanol using the SHAKE (34) and RATTLE (35) algorithms. These bonds were stretched from their equilibrium bond lengths in increments of 0.30 Bohr Radii (0.158 Å) and constrained to those bond lengths while performing AIMD: 1.2 fs time step, 300 K temperature. The mean constraint force () as a function of the constrained bond length (reaction coordinate) was evaluated from a 2-ps AIMD trajectory after equilibration was established (at least a 300-fs period). To determine the free energy of activation, the was integrated over the reaction coordinate using the trapezoidal rule.

Herron et al.

1. Hamnett A (1997) Mechanism and electrocatalysis in the direct methanol fuel cell. Catal Today 38(4):445–457. 2. Arico AS, Srinivasan S, Antonucci V (2001) DMFCs: From fundamental aspects to technology development. Fuel Cells (Weinh) 1(2):133–161. 3. Liu HS, et al. (2006) A review of anode catalysis in the direct methanol fuel cell. J Power Sources 155(2):95–110. 4. Greeley J, Mavrikakis M (2004) Alloy catalysts designed from first principles. Nat Mater 3(11):810–815. 5. Sha Y, Yu TH, Merinov BV, Goddard WA (2012) Prediction of the dependence of the fuel cell oxygen reduction reactions on operating voltage from DFT calculations. J Phys Chem C 116(10):6166–6173. 6. Zhao J, Chan CT, Che JG (2007) Effects of an electric field on a water bilayer on Ag(111). Phys Rev B 75(8):085435. 7. Filhol JS, Neurock M (2006) Elucidation of the electrochemical activation of water over Pd by first principles. Angew Chem Int Ed Engl 45(3):402–406. 8. Jinnouchi R, Anderson AB (2008) Electronic structure calculations of liquid-solid interfaces: Combination of density functional theory and modified Poisson-Boltzmann theory. Phys Rev B 77(24):245417. 9. Skúlason E, et al. (2007) Density functional theory calculations for the hydrogen evolution reaction in an electrochemical double layer on the Pt(111) electrode. Phys Chem Chem Phys 9(25):3241–3250. 10. Janik MJ, Taylor CD, Neurock M (2007) First principles analysis of the electrocatalytic oxidation of methanol and carbon monoxide. Top Catal 46(3-4):306–319. 11. Taylor CD, Wasileski SA, Filhol JS, Neurock M (2006) First principles reaction modeling of the electrochemical interface: Consideration and calculation of a tunable surface potential from atomic and electronic structure. Phys Rev B 73(16):165402. 12. Nørskov JK, et al. (2004) Origin of the overpotential for oxygen reduction at a fuel-cell cathode. J Phys Chem B 108(46):17886–17892. 13. Letchworth-Weaver K, Arias TA (2012) Joint density functional theory of the electrode-electrolyte interface: Application to fixed electrode potentials, interfacial capacitances, and potentials of zero charge. Phys Rev B 86(7):075140. 14. Otani M, Sugino O (2006) First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach. Phys Rev B 73(11):1115407. 15. Ciccotti G, Ferrario M (2004) Blue Moon approach to rare events. Mol Simul 30(11-12): 787–793. 16. Ciccotti G, Kapral R, Vanden-Eijnden E (2005) Blue Moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics. ChemPhysChem 6(9):1809–1814. 17. Mulders T, Kruger P, Swegat W, Schlitter J (1996) Free energy as the potential of mean constraint force. J Chem Phys 104(12):4869–4870.

18. Sprik M, Ciccotti G (1998) Free energy from constrained molecular dynamics. J Chem Phys 109(18):7737–7744. 19. Sprik M (1998) Coordination numbers as reaction coordinates in constrained molecular dynamics. Faraday Discuss 110:437–445. 20. Meijer EJ, Sprik M (1998) Ab initio molecular dynamics study of the reaction of water with formaldehyde in sulfuric acid solution. J Am Chem Soc 120(25):6345–6355. 21. Perdew JP, Burke K, Ernzerhof M (1996) Generalized gradient approximation made simple. Phys Rev Lett 77(18):3865–3868. 22. Karp EM, Silbaugh TL, Crowe MC, Campbell CT (2012) Energetics of adsorbed methanol and methoxy on Pt(111) by microcalorimetry. J Am Chem Soc 134(50): 20388–20395. 23. Greeley J, Mavrikakis M (2002) A first-principles study of methanol decomposition on Pt(111). J Am Chem Soc 124(24):7193–7201. 24. Greeley J, Mavrikakis M (2004) Competitive paths for methanol decomposition on Pt(111). J Am Chem Soc 126(12):3910–3919. 25. Otani M, et al. (2008) Structure of the water/platinum interface—A first principles simulation under bias potential. Phys Chem Chem Phys 10(25):3609–3612. 26. Henkelman G, Uberuaga BP, Jonsson H (2000) A climbing image nudged elastic band method for finding saddle points and minimum energy paths. J Chem Phys 113(22): 9901–9904. 27. Morikawa Y (2001) Adsorption geometries and vibrational modes of C2H2 on the Si(001) surface. Phys Rev B 63(3):033405. 28. Morikawa Y, Ishii H, Seki K (2004) Theoretical study of n-alkane adsorption on metal surfaces. Phys Rev B 69(4):041403. 29. Haynes WH, ed (2011) CRC Handbook of Chemistry and Physics (CRC, New York), 92nd Ed. 30. Troullier N, Martins JL (1991) Efficient pseudopotentials for plane-wave calculations. Phys Rev B Condens Matter 43(3):1993–2006. 31. Vanderbilt D (1990) Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys Rev B Condens Matter 41(11):7892–7895. 32. Monkhorst HJ, Pack JD (1976) Special points for Brillouin-zone integrations. Phys Rev B 13(12):5188–5192. 33. Methfessel M, Paxton AT (1989) High-precision sampling for Brillouin-zone integration in metals. Phys Rev B Condens Matter 40(6):3616–3621. 34. Ryckaert JP, Ciccotti G, Berendsen HJC (1977) Numerical-integration of Cartesian equations of motion of a system with constraints—Molecular-dynamics of N-alkanes. J Comput Phys 23(3):327–341. 35. Andersen HC (1983) Rattle —A velocity version of the shake algorithm for moleculardynamics calculations. J Comput Phys 52(1):24–34.

Herron et al.

PNAS | Published online August 8, 2016 | E4945

PNAS PLUS

Northwest National Laboratory (PNNL); the Center for Nanoscale Materials (CNM) at Argonne National Laboratory; and the National Energy Research Scientific Computing Center (NERSC). EMSL is sponsored by the Department of Energy’s Office of Biological and Environmental Research located at PNNL. CNM and NERSC are supported by the US Department of Energy, Office of Science, under Contracts DE-AC02-06CH11357 and DE-AC02-05CH11231, respectively.

ENGINEERING

ACKNOWLEDGMENTS. Work at University of Wisconson–Madison has been supported by the US Department of Energy - Basic Energy Sciences (DOE-BES), Division of Chemical Sciences, Grant DE-FG02-05ER15731, and was initiated by National Science Foundation East Asia and Pacific Summer Institutes Travel Grant 1014611. The computational work was performed, in part, using supercomputing resources at the following institutions: the Environmental Molecular Sciences Laboratory (EMSL), a national scientific user facility at Pacific