Theoretical Study of Proton Coupled Electron Transfer Reactions: The Eﬀect of Hydrogen Bond Bending Motion Yang Liu, Hao Liu, Kai Song, Yang Xu, and Qiang Shi* Beijing National Laboratory for Molecular Sciences, State Key Laboratory for Structural Chemistry of Unstable and Stable Species, Institute of Chemistry, Chinese Academy of Sciences, Zhongguancun, Beijing 100190, China ABSTRACT: We investigate theoretically the eﬀect of hydrogen bond bending motion on the proton coupled electron transfer (PCET) reaction, using a model system where an intramolecular hydrogen-bonded phenol group is the proton donor. It is shown that, in a two-dimensional (2D) model of the PCET reaction, the bending and stretching vibrational motions are separated, and due to the hydrogen bond conﬁguration and anharmonicity of the potential energy surface, the bending vibration can play a role in the PCET reaction. The results are also compared with two diﬀerent sets of one-dimensional models (1D-linear and 1D-curved). Due to contributions of the bending motion, the rate constants in the 2D model are larger than those in the 1D-linear model, although the diﬀerences between the total rate constants and KIEs for 2D and 1D models are not major. Results from the 1Dcurved model lie between the 2D- and 1D-linear models, indicating that it can include some eﬀect of bending motion in reducing the potential energies along the reaction path. tally.4,12,14,34−45 In these systems, the phenol group, as the proton donor, forms an intramolecular hydrogen bond with the proton acceptor group, which has allowed experimentalists to control the proton donor and acceptor distances by molecular design. Due to the relatively simple structure of these model systems, they have also been intensively investigated in recent theoretical studies.13,43,45−47 For example, in the study of Auer et al.,46 the rate constants of single and double proton transfer were calculated using theoretical expressions of a concerted electrochemical PCET process. In the work by Markle et al.,47 concerted PCET reactions of three similar molecules with diﬀerent proton donor−acceptor distances were compared. In their studies, the proton was treated as a one-dimensional (1D) particle, and the total rate constant including all contributions of R was calculated using 1D models. The above theoretical studies have provided important insights into the eﬀect of hydrogen bond distance ﬂuctuations and the role of excited vibronic states on PCET reactions.46,47 However, these works have neglected the three-dimensional nature of the proton motion by assuming a 1D description of the proton coordinate. As shown in Scheme 1, the bending and stretching motions of a hydrogen bond represent two diﬀerent degrees of freedom. Both of them can aﬀect the overlap of wave functions and then contribute to the PCET reaction, which is related to the angle of the hydrogen bond. If the hydrogen bond angle of D−H−A is close to 180° (Scheme 1a), the vibrational wave function overlap between the vibronic states

1. INTRODUCTION Proton coupled electron transfer (PCET) reactions are ubiquitous in many chemical, biological, and electrochemical processes.1−8 A typical example of PCET can be found in photosystem II (PS II),9−11 where tyrosine-Z is oxidized by P680+ through a long-range electron transfer (ET), together with the phenolic proton transferred to a nearby histidine residue. In this reaction, the electron and proton are likely transferred in a single kinetic step to avoid a high energy stable intermediate state.10 One important problem in the PCET reaction mechanism is related to the eﬀect of hydrogen bond on the coupled motion of the electron and proton.12−14 In this work, we conﬁne our discussion to the case of concerted PCET, where the ET and proton transfer (PT) occur simultaneously. Many studies have shown that, for concerted PCET, ﬂuctuations of the proton donor−acceptor separation RDA (which will be denoted as R later in this work) strongly aﬀect the vibrational wave function overlap between diﬀerent vibronic states of the proton degree of freedom.6,15−20 Because of this, conﬁgurations with small donor−acceptor distance R, although not signiﬁcantly populated, may have a relatively large contribution to the total PCET rate constant.21,22 This eﬀect often leads to a much smaller kinetic isotope eﬀect (KIE) than that calculated using the overlap of vibrational wave functions at the equilibrium proton donor−acceptor distance.23 Systems like PS II are often too complicated for detailed studies of the reaction mechanism and the role of hydrogen bonds in PCET reactions. Thus, many synthetic model systems have been designed to perform such studies.24−33 Recently, molecules with intramolecular hydrogen-bonded phenols have been designed to study the PCET mechanisms experimen© 2015 American Chemical Society

Received: March 27, 2015 Revised: May 21, 2015 Published: May 21, 2015 8104

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

2. THEORY In this work, we use the rate theory for concerted PCET6,22 that is based on Fermi’s golden rule (FGR).53,54 The perturbative treatment using the FGR is based on the assumption of nonadiabatic reaction that is consistent with recent studies of similar systems,46,47 but it may become invalid when the reaction becomes adiabatic. Such types of theories were ﬁrst introduced by Cukier based on nonadiabatic transitions between multiple vibronic states with a quantized proton degree of freedom.2,55 An important generalization of these initial works is to include explicitly the coupling of the proton coordinate to the proton donor−acceptor motion.6,15−20 A two electronic state model is used to describe the concerted PCET.6,22 The potential energy surfaces on the donor and acceptor electronic states are denoted as Vd,a(q, R) as a function of the proton coordinate q and the proton donor−acceptor distance R, where d and a denote the electron donor and acceptor states, respectively. Here, we have retained the multidimensional nature of the proton coordinate q while assuming a 1D description for the proton donor−acceptor motion. The eﬀect of the solvent degrees of freedom on the motion of the electron and proton is treated as a Gaussian process,6,22 which is similar to that in electron transfer reactions.56,57 It is also assumed that there is a separation of time scales between the PCET reaction and the ﬂuctuations of the proton donor− acceptor distance R, such that we can calculate the charge transfer rate at a given proton donor−acceptor distance R and sum over the rate constants with diﬀerent R values.46,47 This static approximation has been tested in solvable models that include the dynamic eﬀects of the R mode, and was found to be reasonably accurate as long as the characteristic frequency of the R mode is small compared to the proton vibrational frequency.23,58 For a ﬁxed donor−acceptor separation R, the PCET rate constant can be calculated as2,6

Scheme 1. Schematic Depictions Showing the Eﬀect of Bending Motion of the Hydrogen Bond

on the donor and acceptor potential energy surfaces (PESs) will be mainly caused by the stretching vibration, and the 1D description is suﬃcient. However, in many cases, the hydrogen bond angle of D−H−A is much smaller than 180°, especially in cases of intramolecular hydrogen bonds. Thus, the hydrogen bond bending vibration may also play an important role in the PCET reaction (Scheme 1b), and a higher dimensional description of the proton degree of freedom is needed. In this work, we investigate the eﬀect of hydrogen bond bending motion on the concerted PCET reactions beyond the 1D models used in recent studies,46,47 by including explicitly the quantized bending degree of freedom. In order to achieve this, the transferring proton should be treated at least as a twodimensional (2D) particle. Regarding the description of the proton motion, we note that there are previous studies of adiabatic proton transfer reactions in the condensed phase that have employed multidimensional descriptions of the proton motion,48,49 and 1D models are also widely used to study the PCET reactions using the molecular dynamics with quantum transitions (MDQT),50 the ring polymer molecular dynamics (RPMD), 51 and the mixed quantum-classical Liouville approach.52 The model system studied in this work is HOArpy (shown in Scheme 2), which has been recently studied Scheme 2. PCET Reaction in the HOAr-Py Molecule

k(R ) =

∑ Pμ,d(R) ∑ kμν(R) μ

(1)

ν

where Pμ,d is the equilibrium population of vibrational state |μ⟩ on the donor PES and kμν is the transition rate constant from state |μ, d⟩ on the donor PES to state |ν, a⟩ on the acceptor PES. Both P and k are considered as functions of the proton donor−acceptor distance R. The rate constant kμν(R) is calculated as2,6 kμν(R ) ≈

experimentally,14,40,44 and theoretically using the 1D model.47 Due to the conjugate structure of the molecule, we can describe the proton degree of freedom with a simpliﬁed 2D model and study the eﬀect of hydrogen bond bending motion. To better understand the eﬀect of hydrogen bond bending motion, calculations are also performed using 1D models for comparison. The remainder of this paper is arranged as follows. The theoretical model that takes into account the 2D proton coordinates and theories used to calculate PCET rate constants are brieﬂy presented in section 2. The computational details and some necessary approximations to obtain the rate constants are presented in section 3. Results using the 2D model are shown in section 4 and compared with those from the 1D models. Finally, conclusions and discussions are made in section 5.

|Vμν(R )|2 ℏ

π λμν(R )kBT

⎡ (ΔG (R ) + λ (R ))2 ⎤ μν μν ⎥ exp⎢ − 4λμν(R )kBT ⎢⎣ ⎥⎦

(2)

In the above eq 2, Vμν(R) is the vibronic coupling between the vibronic states |μ, d⟩ and |ν, a⟩ Vμν(R ) ≈ Vel(R )Sμν(R )

(3)

where Vel(R) is the electronic coupling between the donor and acceptor states and Sμν(R) is the overlap between the proton vibrational wave functions |μ, d⟩ and |ν, a⟩; ΔGμν(R) is the energy diﬀerence between vibronic states |μ, d⟩ and |ν, a⟩ ΔGμν(R ) = ΔG 0(R ) + [Eν , a(R ) − Eμ , d(R )] 8105

(4)

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B where ΔG0(R) is the (adiabatic) energy diﬀerence between the minima of the reactant and product potential energy surfaces, Eμ,d(R) and Eν,a(R) are the energies (relative to the energy minima of the PES) of the vibrational states, and λμν(R) is the reorganization energy associated with the pair of vibronic states |μ, d⟩ and |ν, a⟩. The above eqs 1−4 are used to calculate the PCET rate constants at a ﬁxed R, while the total rate constant is obtained by averaging k(R) over P(R), which is the thermal distribution of R: k=

∫ P(R)k(R) dR

point energy calculations are performed by DFT calculations with the B3LYP functional and 6-31G(d,p) basis set in the gas phase using Gaussian 09.59 The reduced (neutral) and oxidized (cation) structures of the molecules are ﬁrst optimized, and the average structure of them is used as the approximate transition state for the concerted PCET reaction.47 Using the average of the neutral and cation structures as the approximate transition state is a compromise due to the diﬃculty to obtain the real transition state conﬁguration, which has been discussed in ref 47. After that, all of the atoms except for the transferring proton are ﬁxed at the transition state structure, and the position of the proton is optimized. With the ﬁxed proton donor−acceptor distance R, we then use the plane deﬁned by the averaged position of the optimized transferring proton in the neutral and cation states and the positions of the ﬁxed proton donor (O) and acceptor (N) atoms to calculate the 2D potential energy surface. By comparing the structures of the neutral and cation states, the proton moves primarily in this plane, such that treating the proton as a 2D particle is a reasonable approximation for the system considered in this work. The 2D potential energy surfaces of the proton on the neutral and cation states are calculated on a 40 × 40 grid. The 1D potential curves are calculated from the 2D potential energy surfaces using two diﬀerent methods. In the ﬁrst method, we identify the two potential energy minima on the neutral and cation states, and the 1D potential energy curve is calculated along the line connecting the two potential energy minima.46 This method is then called the 1D-linear method. In the second method to obtain the 1D potential energy curve, we use the method that has been adopted in ref 47 where the 1D coordinate for proton motion is deﬁned as r = 1/2(dOH − dNH), where dOH and dNH are distances between the proton and donor/acceptor atoms. The 1D potential energy curve is then obtained by ﬁnding the minimum energy points on 2D potential energy surfaces by changing a series of r values. The potential energy curve deﬁned in this way is termed as the 1Dcurved model. After the 2D and 1D potential energy surfaces at a speciﬁc donor−acceptor distance R are obtained, the vibrational wave functions and energy levels are calculated using the discrete variable representation (DVR) method.60 The calculated parameters (Pμ, Sμν, Eμ/ν) in eqs 2−4 are then used to obtain k(R). To calculate the classical and quantum distributions of the R mode in eqs 7−9, we perform normal mode calculations on the optimized structure of the reactant (the neutral state) to obtain the vibrational frequencies ωk, reduced masses Mk, and the normal mode unit vectors ek. The coeﬃcients ck’s in eq 6 are then calculated as ck = (RA − RD)·ek.46 k(R) is calculated at seven diﬀerent values of R. To perform the numerical integration in eq 5, the k(R) curve is obtained using cubic spline interpolation.61

(5)

There are several ways in the literature to calculate P(R), which include the harmonic approximation based on an ab initio calculation,47 or the method based on normal mode approximation.46 In this work, we apply an extension of the normal mode approximation method proposed by Auer et al.,46 which is based on a decomposition of δR, the deviation of R from its equilibrium value R0, as a linear combination of the normal modes: 3N

δR = R − R 0 ≈

∑ c kQ k

(6)

k=1

The Boltzmann distribution of R is given by a Gaussian function ⎡ (R − R )2 ⎤ 1 0 ⎥ exp⎢ − 2 2 2π ⟨δR ⟩ 2⟨δR ⟩ ⎦ ⎣

P(R ) =

(7)

If the vibrational modes are treated classically, ⟨δR ⟩ can be calculated as46 2

3N

⟨δR2⟩cl =

∑ k=1

ck 2 Mkβωk 2

(8)

where ωk and Mk are the frequency and eﬀective mass of the kth normal mode. If the vibrational modes are to be treated quantum mechanically, we have 3N

⟨δR2⟩qm =

∑ k=1

ℏck 2 2Mkωk tanh(β ℏωk /2)

(9)

3. COMPUTATIONAL DETAILS The parameters Pμ, Sμν, and Eμ/ν in eqs 2−4 (in main text) are functions of the proton donor−acceptor distance R, and are calculated by solving the time-independent Schö dinger equation at a given R. Since the main goal of this work is to investigate the eﬀect of hydrogen bond bending motion on the PCET reaction, the other parameters Vel, λ, and ΔG0 are assumed to be insensitive to the hydrogen bond donor− acceptor distance R,46,47 and are obtained from previous experimental studies.40 More speciﬁcally, for the reaction with the oxidant [Fe(5,5′-Me2bpy)3]3+ in the MeCN solution, the reorganization energy λ is taken to be 23.2 kcal/mol, and ΔG0 ≈ 0 kcal/mol. Vel is assumed to be 36.2 cm−1, in order to reproduce the experimental rate constant40 of 5.8 × 105 s−1 using the 2D model and the quantum distribution of R. The potential energy surfaces on the reactant and product states are calculated using DFT calculations. For simplicity, the tert-butyl groups of the HOAr-py molecule are replaced with methyl groups. All of the structural optimizations and single-

4. RESULTS AND DISCUSSION We ﬁrst calculate the potential energy surfaces and k(R) at diﬀerent donor−acceptor distances. The equilibrium proton donor−acceptor distances R optimized on the neutral and cation states are similar, which are about 2.55 Å. Since k(R) may become much larger at shorter distances, smaller R values should also be considered. As will be shown later (see section 4.3), the case with R = 2.45 Å contributes signiﬁcantly to the total rate constant in both the classical and quantum treatments of the proton donor−acceptor motion. Therefore, in the 8106

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

Figure 1. (a, b) 2D potential energy surface of the neutral (a) and cation (b) states at equilibrium R = 2.55 Å. (c, d) 2D potential energy surface of the neutral (c) and cation (d) states at compressed R = 2.45 Å.

following, we will show the results in detail at the equilibrium distance R = 2.55 Å, and at the compressed distance R = 2.45 Å. Subsequently, the total rate constant and KIE are calculated, and the eﬀect of temperature is investigated. 4.1. Results at the Equilibrium R = 2.55 Å. The 2D potential energy surfaces of the neutral (reduced) and cation (oxidized) states are shown in Figure 1a and b. It can be seen that the potential energy surfaces of the two diﬀerent states tilt toward diﬀerent directions in the neutral and cation states, which is similar to the Dushinsky rotation eﬀect investigated in our previous work.58 Thus, it is expected that the bending motion may also play a similar role as the proton donor− acceptor motion; i.e., excitation of the bending vibrational motion will also increase the PCET reaction rate. However, it should be noted that the vibrational frequency of the bending motion is usually much higher, and the anharmonicity of the PES is not considered in ref 58. As stated in the computational details, we use two diﬀerent methods to deﬁne the 1D potential curves from the 2D potential surface. The 1D-linear proton coordinate is deﬁned along a straight line (blue line in Figure 2) passing through the potential energy minima (denoted as blue dots in Figure 2) of the two potential energy surfaces shown in Figure 1a and b. The 1D-linear coordinate is not parallel to the line between O and N atoms (X-axis in Figure 2), which should be due to the diﬀerences between the OH and NH bond lengths and the diﬀerent hydrogen bond angles in the neutral and cation states. The 1D-curved coordinate is, however, not uniquely deﬁned, since they are slightly diﬀerent when calculated on the 2D neutral and cation PESs (the black and red lines in Figure 2).

Figure 2. Deﬁnition of the reaction coordinates in the two diﬀerent 1D models, on the plane of the 2D model. The two potential minima on the neutral and cation states are shown as blue dots. The reaction coordinate of the 1D-linear model is shown as the blue line. Proton positions that are used to deﬁne the 1D-curved model on the neutral and cation states are shown as black and red curves, respectively.

Strictly speaking, the overlap integral between the neutral and cation states should be calculated using the same set of curved coordinates. However, since they are very close to each other, we have neglected their diﬀerences in our calculations. The calculated potential energy curves using the 1D-linear and 1D-curved coordinates are then shown in Figure 3 (blue curves). The 1D-curved potential energy curve is generally lower than the 1D-linear one, since the 1D-curved coordinate follows the lowest energy path on the 2D PES. For example, 8107

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

the stretching motion of the NH bond, and is mixed with the proton transfer state. The excited vibrational state showing proton transfer characters has been observed and discussed previously in the recent theoretical calculations using 1D models.28,43,46,47 This eﬀect is also shown in the 1D wave functions, as depicted in Figure 3. Here, the ﬁrst excited state on the cation potential energy curve shows signiﬁcant proton transfer character. The wave functions from the 1D-linear and 1D-curved potential energy curves are also very similar. Table 1 lists the proton vibrational energy diﬀerences between the lowest three excited states and ground state in the 2D and 1D models. The vibrational frequencies of OH in the neutral and NH in the cation states are also calculated using Gaussian 0959 with all of the atoms ﬁxed at the approximate transition state except for the transferring proton. The bending and stretching vibrational frequencies for the OH bond (neutral state) are 1483 and 2853 cm−1, and those for the NH bond (cation state) are 1370 and 2771 cm−1, respectively. We note that the 0−1 and 0−2 transition energies on the 2D potential energy surfaces are all much smaller than the normal-mode frequencies, indicating that the anharmonicity is signiﬁcant. The results in Table 1 also show that the energy diﬀerences for the cation state are much smaller than those for the neutral state, which is due to the smaller barrier and stronger anharmonicity of the cation PES. To further investigate the reason, the natural bond orbital (NBO) atomic charges are calculated for each state. In the neutral state, the negative charge is localized on the phenol ring, especially on the oxygen atom (−0.702), and the nitrogen atom is also negatively charged, with slightly smaller charges (−0.526). In the cation state, the molecule becomes positive, but the oxygen atom (−0.534), as the hydrogen bond acceptor, still has more negative charge than the nitrogen atom (−0.483). Thus, it is expected that the NH bond in the cation state is weaker than the OH bond in the neutral state, which results in stronger anharmonicity of the cation PES. Vibrational wave functions of deuteron are also calculated and shown in Figure 5. Comparing with H, the vibrational wave functions of D are more localized and the energy levels of D are smaller than H (results of the energy levels not shown for D).

Figure 3. (a, b) 1D potential energy curves (blue) and the proton vibrational wave function of the ground (black) and the ﬁrst (red) and second (green) excited states in the 1D-linear model on the neutral (a) and cation (b) states. (c, d) The proton vibrational wave function of the ground (black) and the ﬁrst (red) and second (green) excited states in the 1D-curved model on the neutral (c) and cation (d) states. All of the calculations are performed at the equilibrium O−N distance R = 2.55 Å.

the potential energy barriers in the 1D-linear model (10.7 kcal/ mol on the neutral PES and 7.4 kcal/mol on the cation PES) are slightly higher than those in the 1D-curved model (10.0 kcal/mol on the neutral PES and 6.7 kcal/mol on the cation PES). It is thus expected that the latter should give slightly larger vibrational wave function overlaps and reaction rates. The time-independent Schrödinger equation is then solved to obtain the vibrational energy levels and wave functions. Parts a−h of Figure 4 depict the 2D proton vibrational wave functions on the neutral and cation PESs for the lowest four vibrational states. On the neutral PES, the ﬁrst, second, and third excited states correspond to the bending, stretching, and ﬁrst bending overtone states of OH. On the cation PES, the ﬁrst (Figure 4f) excited state corresponds to the bending vibration of NH, while the second excited state (Figure 4g) shows clearly the characteristics of a proton transfer state, which is caused by the anharmonicity of the potential energy surface. The second excited state is also mixed with the bending motion to some extent. The third one (Figure 4h) is related to

Figure 4. 2D proton vibrational wave functions of the ground and ﬁrst three excited states, on the neutral (a−d) and cation (e−h) potential energy surfaces, at the equilibrium O−N distance R = 2.55 Å. 8108

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

Table 1. Proton Vibrational Energy Diﬀerences between the First Three Excited States and the Ground State at the Equilibrium R = 2.55 Åa 2D ΔE (neutral) ΔE (cation) a

1D-linear

1D-curved

0−1

0−2

0−3

0−1

0−2

0−3

0−1

0−2

0−3

1416 1176

2202 1672

2777 2363

2093 1522

2990 2101

4056 3315

1885 1428

2905 2060

3898 3214

The unit of the energy diﬀerences is cm−1.

Figure 5. 2D vibrational wave functions of the ground and ﬁrst three excited states, on the neutral (a−d) and cation (e−h) potential energy surfaces for deuteron transfer, at the equilibrium O−N distance R = 2.55 Å.

Table 2. Reaction Rate Constant, KIE, Vibrational Function Overlaps, and Contributions to the Total Rate Constant, for Proton and Deuteron Coupled Electron Transfer Calculated at the Equilibrium R = 2.55 Åa Sμ,ν2 (percentage) b

k

KIE

(0, 0)

(0, 1)

(0, 2)

H

1.88 × 105

1.97

D

0.96 × 105

H

2.61 × 105

D

1.22 × 105

H

3.07 × 105

D

1.40 × 105

0.002 (20.6%) 0.00004 (0.6%) 0.0044 (28.6%) 0.0001 (1.6%) 0.004 (23.6%) 0.0001 (1.1%)

0.394 (62.2%) 0.018 (9.1%) 0.378 (55.8%) 0.028 (14.5%) 0.142 (36.8%) 0.0033 (3.9%)

0.453 (12.8%) 0.654 (82.2%) 0.47 (10.6%) 0.61 (73.2%) 0.527 (33%) 0.217 (44.5%)

model 1D linear

1D curved

2D

2.13

2.2

(0, 3)

(0, 4)

0.12

0.025

0.199 (3.8%)

0.097 (0.1%)

0.1 (4.6%) 0.187 (1.4%) 0.274 (25.9%)

0.04 (0.2%) 0.213 (15.7%)

a

H and D represent the proton and deuteron, respectively. A blank entry in the table means that the contribution of that speciﬁc pathway is very small. bRate constants are in units of s−1.

total rate constant); thus, we only include the (0, ν) transitions in the following discussions. We note that the total rate constants listed in Table 2 are calculated by including contributions from both ground and excited vibrational states on the reactant PES. The contributions of each pair of vibronic states (μ, ν) are aﬀected by both the vibrational energy diﬀerence and the vibrational wave function overlap. For the proton transfer reaction, the square of the overlap integral between the ground states S0,02 = 0.004 in the 2D model, which is similar to that in the 1D-curved model, is about twice of that in the 1D-linear model. This is because the 2D vibrational wave functions of the ground states are tilted elliptic distributions (Figure 4a and e),

Similar to the situation of proton transfer, deuteron transfer (DT) states are also observed in excited vibrational states on the cation PES. The diﬀerences of vibrational states between H and D are the origin of the KIE, which will be discussed later. Table 2 shows the calculated total rate constant, square of vibrational wave function overlaps, and contributions to the total rate constant from each pair of vibronic states, for proton (H) and deuteron (D) transfer at R = 2.55 Å. Despite that the ﬁrst excited state (bending vibration) on the reactant PES has a lower vibrational energy in the 2D model, the Boltzmann population of the bending excited state at room temperature is still very small (0.0001 for H and 0.0069 for D), which leads to very small contributions from the k1ν terms (less than 3% of the 8109

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

Figure 6. Same as Figure 4, at the compressed O−N distance of R = 2.45 Å.

excluded, the rate constant in the 2D model (1.94 × 105 s−1) is similar to that in the 1D-linear model (1.88 × 105 s−1). Because the potential energies in the 1D-curved model are lower than those in the 1D-linear model, larger wave function overlaps and smaller energy diﬀerences are obtained in the 1D-curved model. Thus, the 1D-curved model can in some aspect include the eﬀect of the bending motion in reducing the potential energies along the reaction path. For deuteron transfer, the more localized vibrational wave functions lead to much smaller values of S0,02 and S0,12 in the 2D model. As a result, the (0, 2) transition to the DT state and the (0, 3) transition to the stretching vibrational state (mixed with some DT character) on the cation PES become important. In the 1D models, the dominant pathway is the (0, 2) transition to the DT state. The KIEs calculated using diﬀerent models are all around 2. We note that, even at R = 2.55 Å, the KIE around 2 is very close to that at R = 2.45 Å (see later), and the experimental value (2.5 ± 0.6).40 Thus, the KIE itself is not only aﬀected by the donor−acceptor distance R. The anharmonicity of the PES, which leads to much delocalized excited vibrational wave functions and smaller energy levels, also plays an important role in determining the KIE. 4.2. Results at the Compressed R = 2.45 Å. In this subsection, we present the results for the compressed proton donor−acceptor distance R = 2.45 Å. In comparison with the molecular structure at R = 2.55 Å, the optimized OH bond length changes from 1.009 to 1.022 Å, in the neutral state, and the NH bond length changes from 1.056 to 1.070 Å in the cation state. The hydrogen bond angle also changes slightly: from 150.1 to 151.6° in the neutral state and from 143.6 to 145.9° in the cation state. The distance between the optimized proton position in the two states also decreases from 0.623 to 0.471 Å. We note that, although the donor−acceptor distance R is shortened by only 0.1 Å, the proton transfer distance is reduced by about 0.15 Å, which is due to the change of bond lengths and angles. This observation is similar to a previous study by Johannissen et al.13 In their work, the proton donor− acceptor distances in two diﬀerent molecules diﬀer about 0.07 Å, but the proton transfer distances are almost identical within 0.01 Å because one molecule has a more linear intramolecular hydrogen bond. Thus, the change of donor−acceptor distance cannot fully determine the change of proton tunneling distance

which will lead to a larger overlap between two ground states below the line connecting the two potential minima. The diﬀerence of S0,02 results in a smaller contribution of the (0, 0) pathway in the 1D-linear model. Diﬀerences between the 2D and 1D models become obvious when the transitions to the excited vibrational states on the product PES are included. The (0, 1) transition in the 2D model represents the pathway from the ground state on the neutral PES to the bending vibrational state on the cation PES. Because the bending vibration of NH on the cation PES tilts toward the proton equilibrium position on the neutral PES, the (0, 1) overlap is much larger than the (0, 0) one, with S0,12 = 0.142. The (0, 1) transition also makes the largest contribution (36.8%) to the total rate constant. The (0, 2) transition has the largest wave function overlap S0,22 = 0.527 due to transition to the PT state on the cation PES, but its contribution (33%) to the total rate constant is slightly smaller because of the larger vibrational energy diﬀerence. Similarly, although S0,32 = 0.187 is slightly larger than S0,12, contribution of the (0, 3) pathway (1.4%) is very small due to the even larger vibrational energy diﬀerence. Thus, the major contribution of the excited vibrational states (on the cation PES) is attributed to the bending state (due to its low vibrational energy) and the proton transfer state (due to its large wave function overlap). Anharmonicity of the cation PES also plays an important role for the large S0,12 value. As shown in Figure 1b and Figure 4, due to the larger anharmonicity on the cation state, the N−H bending vibrational state has a larger spread toward the potential minimum on the donor PES, which results in larger overlap with the ground vibrational state on the neutral state. In contrast, the neutral PES is less anharmonic, and S1,02 = 0.012 is much smaller than S0,12 = 0.142. Thus, contributions of the bending vibrational state could vary in systems with diﬀerent hydrogen bond conﬁgurations, and may be altered by molecular design.44 Since the bending motion is neglected in the 1D models, the calculated rate constants in the 1D models are smaller than that in the 2D model. In the 1D models, the dominant (0, 1) pathway corresponds to transition from the ground vibrational state on the neutral PES to the proton transfer state on the cation PES, which is similar to the (0, 2) pathway in the 2D model. If the contribution of the bending vibrational state is 8110

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B Table 3. Same as Table 1, at the Compressed R = 2.45 Å 2D ΔE (neutral) ΔE (cation)

1D-linear

1D-curved

0−1

0−2

0−3

0−1

0−2

0−3

0−1

0−2

0−3

1408 1018

1897 1727

2792 2361

1658 1160

2857 2328

4366 3845

1588 1160

2792 2288

4169 3676

due to the variations of bond lengths and hydrogen bond angles. The 2D potential energy surfaces at R = 2.45 Å are shown in Figure 1c and d. As compared with those at R = 2.55 Å, there are no double well features in either the neutral or cation states, and the potential wells also become signiﬁcantly shallower. The 2D vibrational wave functions are shown in Figure 6, and vibrational energies are listed in Table 3. The vibrational wave functions are similar to those in the case of R = 2.55 Å in Figure 4, except for the third excited state on the cation PES, which shows strong mixing of the bending overtone and the proton transfer states. A clearly deﬁned stretching excited state like that in Figure 4h is also not observed in the R = 2.45 Å case. The 1D potential energy curves and the proton vibrational wave functions are provided in Figure 7. Since the energy curves become shallower, the vibrational wave functions also become more delocalized.

The changes of vibrational wave functions and energy levels further aﬀect the wave function overlaps and rate constants, as shown in Table 4. For proton transfer, reducing R leads to a much larger overlap between ground vibrational states (with S0,02 ≈ 0.1), such that the (0, 0) transition dominates the total reaction constant for R = 2.45 Å. The (0, 1) pathway involving the bending mode still contributes 21.3% to the total rate constant, while contributions from the other vibrational excited states become very small. Since contributions from excited vibrational states become less important, the diﬀerences of rate constants between the three models become smaller. For instance, k2D/k1D‑linear decreases from 1.63 to 1.18 when compared to the case of R = 2.55 Å, while k1D‑curved still lies between k1D‑linear and k2D. For deuteron transfer, both S0,02 and S0,12 increase signiﬁcantly, but S0,02 is still very small, and the (0, 1) transition is the dominant pathway (38.1%). The ﬁrst excited vibrational state of the 2D model is related to the bending mode, which is not considered in the 1D models. This leads to a larger k2D/ k1D‑linear value for DT than PT. The KIEs calculated at R = 2.45 Å for the 2D model are almost the same as that at R = 2.55 Å, while an obvious increase of the KIE of about 0.5 is found for the 1D-linear model. Again, this is due to the fact that the 1D models do not include the transition to the ﬁrst bending excited vibrational state, which is the dominant pathway in the 2D model in deuteron transfer. 4.3. The Total Rate Constant. With further decrease of the proton donor−acceptor distance R, the rate constants increase for both the proton and deuteron transfer. However, the shorter R values contribute little to the total rate constant due to their small thermal populations. The calculated rate constant k(R) in the 2D model, the classical and quantum distributions Pcl(R) and Pqm(R), as well as their product P(R) k(R) are shown in Figure 8. The solid black line is a spline interpolation of the scattered k(R) data calculated at seven diﬀerent R values. It can be seen that the quantum distribution is wider than the classical one, while, for either distribution, the

Figure 7. Same as Figure 3, at the compressed O−N distance of R = 2.45 Å.

Table 4. Same as Table 2, at the Compressed R = 2.45 Åa Sμ,ν2 (percentage) b

k

KIE

(0, 0)

(0, 1)

(0, 2)

(1, 0)

H

2.15 × 106

2.46

D

0.87 × 106

H

2.29 × 106

D

1.0 × 106

H

2.55 × 106

D

1.16 × 106

0.09 (73.2%) 0.012 (24.1%) 0.101 (75.1%) 0.017 (29.1%) 0.108 (72.5%) 0.017 (26.2%)

0.57 (21.8%) 0.304 (51%) 0.552 (19.8%) 0.302 (46.4%) 0.443 (21.3%) 0.187 (38.1%)

0.275 (0.3%) 0.456 (11.7%) 0.286 (0.3%) 0.455 (10.2%) 0.245 (1.6%) 0.338 (18.6%)

0.5 (4.6%) 0.157 (8.4%) 0.432 (4.6%) 0.152 (9.4%) 0.196 (3.1%) 0.040 (4.0%)

model 1D linear

1D curved

2D

a

2.29

2.27

H and D represent the proton and deuteron, respectively. bRate constants are in the units of s−1. 8111

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

that the three diﬀerent models give essentially the same temperature dependence of the total rate constants. For example, k2D/k1D‑linear is 1.33 at 275 K and 1.31 at 315 K. This could be attributed to two competitive factors: the proton donor−acceptor R tends to sample shorter values as the temperature increases, which reduces the diﬀerence of rate constants between the 2D and 1D models, and contributions from excited vibrational states become larger at high temperature, which increases the diﬀerence between them.

5. CONCLUSIONS AND DISCUSSIONS In this work, we have considered the eﬀect of hydrogen bond bending motion on a concerted PCET reaction in a model of HOAr-py oxidation. 2D potential energy surfaces at the approximate transition state conﬁgurations are calculated, and the PCET rate theory based on the FGR is applied to obtain the rate constants, KIE, and contributions from each transition pathway. In the 2D model, the bending and stretching vibrational states are found to be separated. The vibrational excited states on the reactant (neutral) PES have little contribution to the total rate constant because of their small thermal populations. However, due to the anharmonicity of the PES and the small hydrogen bond angle, the bending excited state and proton transfer state are found to play important roles in the PCET reaction. For example, for proton transfer, the major contributions at the equilibrium distance R = 2.55 Å come from the (0, 1) transition to the bending state (due to its low vibrational energy) and the (0, 2) transition to the proton transfer state (due to its large wave function overlap). As R is compressed to the small values (e.g., R = 2.45 Å), the (0, 0) transition becomes the dominant pathway, but due to the smaller energy level, the bending vibrational state still plays a role in the PCET reaction. For comparison, we have also performed calculations using two diﬀerent 1D potential models constructed from the 2D PES. The smaller rate constant in the 1D-linear model can further conﬁrm the contribution of the hydrogen bond bending motion. Since the 1D-curved model includes some eﬀect of bending motion, the results calculated in the 1D-curved model are between the 2D- and 1D-linear models. The temperature almost does not change the relative value of rate constants calculated with the three models. The KIEs in all three models are in good agreement with the experimental value. In summary, the 2D study provides more detailed information on the proton motion, including the anharmonic 2D potential energy surfaces and the 2D proton vibrational wave functions. Although the diﬀerences between the total rate constants and KIEs for the 2D and 1D models are not major, we do observe cases where the transition to bending vibrational states makes the dominant contribution to the total rate constant (e.g., R = 2.55 Å for PT and R = 2.45 Å for DT). To further probe this eﬀect, the current work will be extended to other PCET model systems with diﬀerent hydrogen bond conﬁgurations. As we have employed the perturbative FGR to treat the nonadiabtic PCET reactions, several other approaches including the MDQT,50 the RPMD,51 and the mixed quantumclassical Liouville equation52 can also be used to study the eﬀect of hydrogen bond bending motion, and can go beyond the nonadiabatic limit. The current study has also focused on concerted PCET reactions, but based on the calculated potential energy surfaces and wave function overlaps, the hydrogen bond bending motion may also play a role in

Figure 8. Total reaction constant k(R) at diﬀerent R values (black dots) and their spline interpolation (black line), the classical (red line) and quantum (green line) distributions P(R), and their products k(R) P(R) (classical, blue line; quantum, violet line). The results are calculated using the 2D model at 295 K.

range of R between 2.4 and 2.6 Å dominates the contribution to the total rate constant. The total rate constants for the three diﬀerent models are then calculated using eq 5 and are listed in Table 5 for both H Table 5. Total Rate Constant Calculated for the Three Models with Classical and Quantum Mechanical Distribution Functions of Ra 2D kcl(H) kqm(H) kcl(D) kqm(D) a

4.82 5.8 2.10 2.54

× × × ×

1D-linear 105 105 105 105

3.51 4.39 1.54 1.91

× × × ×

105 105 105 105

1D-curved 4.19 5.11 1.83 2.25

× × × ×

105 105 105 105

The temperature is 295 K. The unit of the rate constants is s−1.

and D transfers. The larger rate constants calculated with the quantum mechanical distribution are attributed to the wider distribution of Pqm(R). The total rate constants in the three models and their ratios are between the results calculated at R = 2.55 and 2.45 Å. All of the KIEs calculated using diﬀerent models are around 2.3, which is in good agreement with the experimental value (2.5 ± 0.6).40 The temperature dependence of the total rate constants is calculated using Pqm(R) and provided in Figure 9. It is found

Figure 9. Temperature dependence of total rate constant calculated with the quantum mechanical distribution Pqm(R) for the 1D-linear, 1D-curved, and 2D models. 8112

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B

(16) Cukier, R. I. Theory and Simulation of Proton-Coupled Electron Transfer, Hydrogen-Atom Transfer, and Proton Translocation in Proteins. Biochim. Biophys. Acta 2004, 1655, 37−44. (17) Soudackov, A.; Hatcher, E.; Hammes-Schiffer, S. Quantum and Dynamical Effects of Proton Donor-Acceptor Vibrational Motion in Nonadiabatic Proton-Coupled Electron Transfer Reactions. J. Chem. Phys. 2005, 122, 014505. (18) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. Comparison of Dynamical Aspects of Nonadiabatic Electron, Proton, and ProtonCoupled Electron Transfer Reactions. Chem. Phys. 2005, 319, 93−100. (19) Klinman, J. P. The Role of Tunneling in Enzyme Catalysis of CH Activation. Biochim. Biophys. Acta 2006, 1757, 981−987. (20) Presse, S.; Silbey, R. J. Anomalous Temperature-Isotope Dependence in Proton-Coupled Electron Transfer. J. Chem. Phys. 2006, 124, 124504. (21) Hatcher, E.; Soudackov, A. V.; Hammes-Schiffer, S. ProtonCoupled Electron Transfer in Soybean Lipoxygenase. J. Am. Chem. Soc. 2004, 126, 5763−5775. (22) Hammes-Schiffer, S.; Stuchebrukhov, A. A. Theory of Coupled Electron and Proton Transfer Reactions. Chem. Rev. 2010, 110, 6939− 6960. (23) Edwards, S. J.; Soudackov, A. V.; Hammes-Schiffer, S. Analysis of Kinetic Isotope Effects for Proton-Coupled Electron Transfer Reactions. J. Phys. Chem. A 2009, 113, 2117−2126. (24) Biczók, L.; Gupta, N.; Linschitz, H. Coupled Electron-Proton Transfer in Interactions of Triplet C60 with Hydrogen-Bonded Phenols: Effects of Solvation, Deuteration, and Redox Potentials. J. Am. Chem. Soc. 1997, 119, 12601−12609. (25) Lachaud, F.; Quaranta, A.; Pellegrin, Y.; Dorlet, P.; Charlot, M.F.; Un, S.; Leibl, W.; Aukauloo, A. A Biomimetic Model of the Electron Transfer between P680 and the TyrZ-His190 Pair of PSII. Angew. Chem., Int. Ed. 2005, 44, 1536−1540. (26) Costentin, C.; Robert, M.; Savéant, J.-M. Concerted ProtonElectron Transfers in the Oxidation of Phenols. Phys. Chem. Chem. Phys. 2010, 12, 11179−11190. (27) Bonin, J.; Costentin, C.; Louault, C.; Robert, M.; Routier, M.; Savéant, J.-M. Intrinsic Reactivity and Driving Force Dependence in Concerted Proton-Electron Transfers to Water Illustrated by Phenol Oxidation. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 3367−3372. (28) Westlake, B. C.; Brennaman, M. K.; Concepcion, J. J.; Paul, J. J.; Bettis, S. E.; Hampton, S. D.; Miller, S. A.; Lebedeva, N. V.; Forbes, M. D.; Moran, A. M.; Meyer, T. J.; Papanikolas, J. M. Concerted ElectronProton Transfer in the Optical Excitation of Hydrogen-Bonded Dyes. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 8554−8558. (29) Sjödin, M.; Styring, S.; Wolpher, H.; Xu, Y.; Sun, L.; Hammarström, L. Switching the Redox Mechanism: Models for Proton-Coupled Electron Transfer from Tyrosine and Tryptophan. J. Am. Chem. Soc. 2005, 127, 3855−3863. (30) Magnuson, A.; Berglund, H.; Korall, P.; Hammarström, L.; Åkermark, B.; Styring, S.; Sun, L. Mimicking Electron Transfer Reactions in Photosystem II: Synthesis and Photochemical Characterization of a Ruthenium (II) Tris (bipyridyl) Complex with a Covalently Linked Tyrosine. J. Am. Chem. Soc. 1997, 119, 10720− 10725. (31) Sjödin, M.; Styring, S.; Åkermark, B.; Sun, L.; Hammarström, L. Proton-Coupled Electron Transfer from Tyrosine in a TyrosineRuthenium-Tris-Bipyridine Complex: Comparison with TyrosineZ Oxidation in Photosystem II. J. Am. Chem. Soc. 2000, 122, 3932−3936. (32) Zhang, M.-T.; Hammarström, L. Proton-Coupled Electron Transfer from Tryptophan: A Concerted Mechanism with Water as Proton Acceptor. J. Am. Chem. Soc. 2011, 133, 8806−8809. (33) Mayer, J. M.; Rhile, I. J.; Larsen, F. B.; Mader, E. A.; Markle, T. F.; DiPasquale, A. G. Models for Proton-Coupled Electron Transfer in Photosystem II. Photosynth. Res. 2006, 87, 3−20. (34) Maki, T.; Araki, Y.; Ishida, Y.; Onomura, O.; Matsumura, Y. Construction of Persistent Phenoxyl Radical with Intramolecular Hydrogen Bonding. J. Am. Chem. Soc. 2001, 123, 3371−3372. (35) Irebo, T.; Johansson, O.; Hammarström, L. The Rate Ladder of Proton-Coupled Tyrosine Oxidation in Water: A Systematic Depend-

sequential PCET reactions, which may be subject to further investigations.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +86 10 82616163. Fax: +86 10 82616163. Notes

The authors declare no competing ﬁnancial interest.

■

ACKNOWLEDGMENTS This work is supported by NSFC (Grant No. 21290194), the 973 program (Grant Nos. 2011CB808502 and 2013CB933501), and the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB12020300).

■

REFERENCES

(1) Tommos, C.; Babcock, G. T. Oxygen Production in Nature: A Light-Driven Metalloradical Enzyme Process. Acc. Chem. Res. 1998, 31, 18−25. (2) Cukier, R. I.; Nocera, D. G. Proton-Coupled Electron Transfer. Annu. Rev. Phys. Chem. 1998, 49, 337−369. (3) Mayer, J. M. Proton-Coupled Electron Transfer: A Reaction Chemist’s View. Annu. Rev. Phys. Chem. 2004, 55, 363. (4) Costentin, C.; Robert, M.; Savéant, J.-M. Electrochemical and Homogeneous Proton-Coupled Electron Transfers: Concerted Pathways in the One-Electron Oxidation of a Phenol Coupled with an Intramolecular Amine-Driven Proton Transfer. J. Am. Chem. Soc. 2006, 128, 4552−4553. (5) Huynh, M. H. V.; Meyer, T. J. Proton-Coupled Electron Transfer. Chem. Rev. 2007, 107, 5004−5064. (6) Hammes-Schiffer, S.; Soudackov, A. V. Proton-Coupled Electron Transfer in Solution, Proteins, and Electrochemistry. J. Phys. Chem. B 2008, 112, 14108−14123. (7) Reece, S. Y.; Nocera, D. G. Proton-Coupled Electron Transfer in Biology: Results from Synergistic Studies in Natural and Model Systems. Annu. Rev. Biochem. 2009, 78, 673−699. (8) Migliore, A.; Polizzi, N. F.; Therien, M. J.; Beratan, D. N. Biochemistry and Theory of Proton-Coupled Electron Transfer. Chem. Rev. 2014, 114, 3381−465. (9) Hoganson, C. W.; Babcock, G. T. A Metalloradical Mechanism for the Generation of Oxygen from Water in Photosynthesis. Science 1997, 277, 1953−1956. (10) Meyer, T. J.; Huynh, M. H. V.; Thorp, H. H. The Possible Role of Proton-Coupled Electron Transfer (PCET) in Water Oxidation by Photosystem II. Angew. Chem., Int. Ed. 2007, 46, 5284−5304. (11) Barry, B. A. Proton Coupled Electron Transfer and Redox Active Tyrosines in Photosystem II. J. Photochem. Photobiol., B 2011, 104, 60−71. (12) Sjödin, M.; Irebo, T.; Utas, J. E.; Lind, J.; Merényi, G.; Åkermark, B.; Hammarström, L. Kinetic Effects of Hydrogen Bonds on Proton-Coupled Electron Transfer from Phenols. J. Am. Chem. Soc. 2006, 128, 13076−13083. (13) Johannissen, L. O.; Irebo, T.; Sjödin, M.; Johansson, O.; Hammarström, L. The Kinetic Effect of Internal Hydrogen Bonds on Proton-Coupled Electron Transfer from Phenols: A Theoretical Analysis with Modeling of Experimental Data. J. Phys. Chem. B 2009, 113, 16214−16225. (14) Markle, T. F.; Mayer, J. M. Concerted Proton-Electron Transfer in Pyridylphenols: The Importance of the Hydrogen Bond. Angew. Chem., Int. Ed. 2008, 47, 738−740. (15) Knapp, M. J.; Rickert, K.; Klinman, J. P. TemperatureDependent Isotope Effects in Soybean Lipoxygenase-1: Correlating Hydrogen Tunneling with Protein Dynamics. J. Am. Chem. Soc. 2002, 124, 3865−3874. 8113

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114

Article

The Journal of Physical Chemistry B ence on Hydrogen Bonds and Protonation State. J. Am. Chem. Soc. 2008, 130, 9194−9195. (36) Zhang, M.-T.; Irebo, T.; Johansson, O.; Hammarström, L. Proton-Coupled Electron Transfer from Tyrosine: A Strong Rate Dependence on Intramolecular Proton Transfer Distance. J. Am. Chem. Soc. 2011, 133, 13224−13227. (37) Costentin, C.; Robert, M.; Savéant, J.-M. Adiabatic and NonAdiabatic Concerted Proton-Electron Transfers. Temperature Effects in the Oxidation of Intramolecularly Hydrogen-Bonded Phenols. J. Am. Chem. Soc. 2007, 129, 9953−9963. (38) Rhile, I. J.; Mayer, J. M. One-Electron Oxidation of a HydrogenBonded Phenol Occurs by Concerted Proton-Coupled Electron Transfer. J. Am. Chem. Soc. 2004, 126, 12718−12719. (39) Rhile, I. J.; Mayer, J. M. Comments on How Single and Bifurcated Hydrogen Bonds Influence Proton-Migration Rate Constants, Redox, and Electronic Properties of Phenoxyl Radicals. Angew. Chem., Int. Ed. 2005, 44, 1598−1599. (40) Rhile, I. J.; Markle, T. F.; Nagao, H.; DiPasquale, A. G.; Lam, O. P.; Lockwood, M. A.; Rotter, K.; Mayer, J. M. Concerted ProtonElectron Transfer in the Oxidation of Hydrogen-Bonded Phenols. J. Am. Chem. Soc. 2006, 128, 6075−6088. (41) Benisvy, L.; Bittl, R.; Bothe, E.; Garner, C. D.; McMaster, J.; Ross, S.; Teutloff, C.; Neese, F. Phenoxyl Radicals Hydrogen-Bonded to Imidazolium: Analogues of Tyrosyl D. of Photosystem II: HighField EPR and DFT Studies. Angew. Chem., Int. Ed. 2005, 44, 5314− 5317. (42) Markle, T. F.; Rhile, I. J.; DiPasquale, A. G.; Mayer, J. M. Probing Concerted Proton-Electron Transfer in Phenol-Imidazoles. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 8185−8190. (43) Markle, T. F.; Rhile, I. J.; Mayer, J. M. Kinetic Effects of Increased Proton Transfer Distance on Proton-Coupled Oxidations of Phenol-Amines. J. Am. Chem. Soc. 2011, 133, 17341−17352. (44) Markle, T. F.; Tronic, T. A.; DiPasquale, A. G.; Kaminsky, W.; Mayer, J. M. Effect of Basic Site Substituents on Concerted ProtonElectron Transfer in Hydrogen-Bonded Pyridyl-Phenols. J. Phys. Chem. A 2012, 116, 12249−12259. (45) Schrauben, J. N.; Cattaneo, M.; Day, T. C.; Tenderholt, A. L.; Mayer, J. M. Multiple-Site Concerted Proton-Electron Transfer Reactions of Hydrogen-Bonded Phenols Are Nonadiabatic and Well Described by Semiclassical Marcus Theory. J. Am. Chem. Soc. 2012, 134, 16635−16645. (46) Auer, B.; Fernandez, L. E.; Hammes-Schiffer, S. Theoretical Analysis of Proton Relays in Electrochemical Proton-Coupled Electron Transfer. J. Am. Chem. Soc. 2011, 133, 8282−8292. (47) Markle, T. F.; Tenderholt, A. L.; Mayer, J. M. Probing Quantum and Dynamic Effects in Concerted Proton−Electron Transfer Reactions of Phenol−Base Compounds. J. Phys. Chem. B 2011, 116, 571−584. (48) Laria, D.; Ciccotti, G.; Ferrario, M.; Kapral, R. MolecularDynamics Study of Adiabatic Proton-Transfer Reactions in Solution. J. Chem. Phys. 1992, 97, 378−388. (49) Agarwal, P. K.; Billeter, S. R.; Hammes-Schiffer, S. Nuclear Quantum Effects and Enzyme Dynamics in Dihydrofolate Reductase Catalysis. J. Phys. Chem. B 2002, 106, 3283−3293. (50) Fang, J. Y.; Hammes-Schiffer, S. Proton-Coupled Electron Transfer Reactions in Solution: Molecular Dynamics with Quantum Transitions for Model Systems. J. Chem. Phys. 1997, 106, 8442−8454. (51) Kretchmer, J. S.; Miller, T. F., III. Direct simulation of protoncoupled electron transfer across multiple regimes. J. Chem. Phys. 2013, 138, 134109. (52) Shakib, F.; Hanna, G. An Analysis of Model Proton-Coupled Electron Transfer Reactions via the Mixed Quantum-Classical Liouville Approach. J. Chem. Phys. 2014, 141, 044122. (53) Weiss, U. Quantum dissipative systems, 2nd ed.; World Scientiﬁc: London, 1999. (54) May, V.; Kühn, O. Charge and Energy Transfer Dynamics in Molecular Systems, 3rd ed.; Wiley-VCH: Weinheim, Germany, 2011. (55) Cukier, R. I. Mechanism for Proton-Coupled Electron-Transfer Reactions. J. Phys. Chem. 1994, 98, 2377−2381.

(56) Marcus, R. A. On the Theory of Oxidation-Reduction Involving Electron Transfer. J. Chem. Phys. 1956, 24, 966−978. (57) Marcus, R. A. Electron Transfer Reactions in Chemistry. Theory and Experiment. Rev. Mod. Phys. 1993, 65, 599−610. (58) Zheng, R. H.; Jing, Y. Y.; Chen, L. P.; Shi, Q. Theory of Proton Coupled Electron Transfer Reactions: Assessing the Born-Oppenheimer Approximation for the Proton Motion Using an Analytically Solvable Model. Chem. Phys. 2011, 379, 39−45. (59) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (60) Colbert, D. T.; Miller, W. H. A Novel Discrete Variable Representation (DVR) for Quantum Mechanical Reactive Scattering via the S-Matrix Kohn Method. J. Chem. Phys. 1992, 96, 1982−1991. (61) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992.

8114

DOI: 10.1021/acs.jpcb.5b02927 J. Phys. Chem. B 2015, 119, 8104−8114