Transition Flux Formula for the Electronic Coupling Matrix Element

Mar 31, 2015 - Here we present a new derivation which is based on the Golden Rule approach. ... The tunneling flux theorem(1) connects the transition ...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Transition Flux Formula for the Electronic Coupling Matrix Element Muhammad A. Hagras and Alexei A. Stuchebrukhov* Department of Chemistry, University of California, One Shields Avenue, Davis, California 95616, United States S Supporting Information *

ABSTRACT: The transition flux formula for the coupling matrix element of longdistance electron transfer reactions is discussed. Here we present a new derivation which is based on the Golden Rule approach. The electronic Franck−Condon factor that appears in the multielectronic formulation of the coupling element is discussed using the concept of tunneling time. An application of the tunneling flux theory to electron transfer reactions in a model system based on the low-potential heme and high-potential heme (heme bL)/(heme bH) redox pair of ubiquinol:cytochrome c oxidoreductase complex is described; the results are compared to those obtained by measuring energy splitting of the donor/acceptor multielectronic states and the direct calculation method.

1. INTRODUCTION The tunneling flux theorem1 connects the transition flux of an electron transfer reaction to the coupling matrix element. The connection is as follows: VDA = −ℏ

∫∂Ω

(d s ⃗ ·J ⃗ )

VDA = −ℏ

where VDA is the coupling matrix element, J ⃗ is the electronic transition flux, and ∂ΩD is the dividing surface between the donor and acceptor complexes, located in regions ΩD and ΩA, respectively. In addition to the coupling matrix element, the spatial distribution of flux J(⃗ r)⃗ carries information about how the charge transition actually occurs in space, i.e., describes tunneling pathways. The tunneling calculations involve two (nonorthogonal) diabatic electronic states |D⟩ and |A⟩, corresponding to localized charge on the donor and on the acceptor respectively, and the transition flux J(⃗ r)⃗ is the matrix element between states A and D of the current density operator:

J ⃗ ( r ⃗) = −i⟨A| j ⃗ ̂ ( r ⃗)|D⟩

kET =

(3)

2 ⎡ (ΔG + λ)2 ⎤ 2π ⟨VDA ⟩ 0 ⎥ exp⎢ − 4λkBT ⎦ ℏ 4πλkBT ⎣

(4)

where λ is the reorganization energy, ΔG0 is the free energy change (−ΔG0 is the driving force), VDA is the coupling matrix element discussed above, and the brackets stand for the statistical averaging over possible variations in the transition state configuration of the system, i.e., the configuration at which the resonance between the donor and acceptor diabatic states occurs; the variations are due to the different realizations of the transition state that occur along the dynamic trajectory of the system.1b One of the most recent applications of the above approach to electron transfer between FeS clusters in a protein (complex I) of the respiratory chain is described in ref 3.

(2)

The formulation is the same for one- or multielectronic description of the system. In order to describe the tunneling pathways in a more convenient language of chemical structures, one can introduce the so-called interatomic currents Jab, which are just coarse-grained currents flowing between atoms. The total current through an atom is proportional to the probability that the tunneling electron will pass through this atom during the tunneling jump; thus, the concept allows one to determine the atoms of the molecular structure (such as a protein) that are involved in the tunneling transition. Moreover, interatomic currents Jab are also connected to the coupling matrix element in a way similar to that of spatial flux J,⃗ namely, © 2015 American Chemical Society

Jab

where the sum of interatomic currents crossing the dividing surface is calculated, instead of the surface integral in eq 1. In practice, it is the latter formulation that is preferred over that involving more detailed flux J(⃗ r)⃗ , because the interatomic currents can be calculated directly from the Hamiltonian of the system, bypassing evaluation of J(⃗ r)⃗ altogether. Both interatomic currents Jab and current density J(⃗ r)⃗ provide full information about the tunneling process, and eventually relate this picture to the rate of electron transfer reaction by the Marcus−Levich−Dogonadze equation2

(1)

D

∑ a ∈ΩD, b ∈ΩA

Special Issue: John R. Miller and Marshall D. Newton Festschrift Received: December 19, 2014 Revised: March 30, 2015 Published: March 31, 2015 7712

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B Previously,1 the tunneling flux formulation was given by using the time-dependent picture of an oscillating wave function at a “frozen” transition state of the systeman artificial construct in which one assumes that the nuclear configuration of the system is fixed at a resonant state, and the total electronic wave function is oscillating between the donor and acceptor states with frequency VDA/ℏ ⎛V ⎞ ⎛V ⎞ |Ψ(t )⟩ = cos⎜ DA t ⎟|D⟩ − i sin⎜ DA t ⎟|A⟩ ⎝ ℏ ⎠ ⎝ ℏ ⎠

Franck−Condon density of vibronic states sometimes used to derive eq 4. As such, no nuclear reorganization energy is contained in VDA; this is purely electronic coupling.) The quantum flux in the above equation is4 j⃗ =

(5)

Ψ(t ) = φD +

j ⃗ ·ds ⃗

j⃗ =

2π VDA 2ρ = ℏ

J⃗ =

D

j ⃗ ·ds ⃗

E

E

ℏ (φ ∇φ − φA∇φD) 2m D A

(12)

We assume that the spatial parts of the donor and acceptor wave functions are real-valued. For time-dependent coefficients CE(t), the perturbation theory gives the following: t 1 VDA e−iE(t − t ′)/ ℏ dt ′ iℏ 0 V (cos(Et /ℏ) − 1) − i sin(Et /ℏ) = DA E/ℏ ℏ



CE(t ) =

(13)

So that, Im CE = −

VDA sin(Et /ℏ) E/ℏ ℏ

(14)

Replacing the sum over continuum energies E in eq 11 with an integral, we have

∑ Im CE = − E

VDA ℏ

∫ ρ dE

sin(Et /ℏ) = −πVDAρ (E /ℏ)

(15)

We now recall the expression for the rate in eq 8, and using eqs 11 and 15, we find 2π VD A 2ρ = −(2πVDAρ) ℏ

∫∂Ω

D

J ⃗ ·ds ⃗

(16)

Cancellation of similar terms on both sides gives the connection:

VDA = −ℏ

(7)

∫∂Ω

∑ Im CE(φD∇φA − φA∇φD) = J ⃗ ·2 ∑ Im CE

where

∫∂Ω

D

The rate of decay is given by the Golden Rule; therefore, at short times, kETt ≪ 1, we have kET =

ℏ m

(11)

where ΩD is the region in space that encloses the donor complex, and ∂ΩD is its boundary. The decay is exponential, so that PḊ = −kETe−kETt

(10)

The summation is over the continuum of the acceptor states, with energies E counted from that of the donor state (i.e., the energy of the donor state is taken to be zero). For the wave function of such a form, we have

(6)

D

∑ CE(t )φEA E

2. TUNNELING FLUX AND THE GOLDEN RULE A tunneling transition between a localized donor state and a continuum of acceptor states can be formally described as follows. (By a “continuum” of acceptor states, we actually mean only one purely electronic state smeared in energy due to nuclear motions of the system. No vibrational Franck−Condon effects are present here, since we consider nuclear motion classically. All nuclear reorganization is entering separately into the exponential factor in eq 4.) The decay of the donor state is related to quantum flux by the continuity equation

∫∂Ω

(9)

and at short times kETt ≪1, when the population of the donor state has not yet changed significantly from its initial value (i.e., unity), the time-dependent wave function has the form

The analysis of periodic changes of populations of the donor and acceptor states occurring during the time evolution of the system results in the above formulation connecting transition flux to VDA. Such a picture of oscillating wave function of course never occurs in reality; the real process is best described in terms of a sequence of Landau−Zener transitions,1b,2 or in the framework of the equivalent Golden Rule that gives directly the rate of the reaction, and exponential decay of the initial (donor) state. All formulations, however, contain the same coupling element VDA, and of course, it should not matter how one calculates it. However, it is more convenient, if not more convincing, to formulate theory in such terms that are more adequate to real physics of the system. Here we present a new derivation which is based on the Golden Rule framework, instead of eq 5. In addition, some discussion will be added to that existing in the literature regarding the multielectronic nature of the transition coupling matrix element, and the electronic Franck−Condon factor that appears in the expression for the coupling when all electrons of the system are considered. The appearance of such a factor was first described by Marshall Newton and then resurfaced in our work on the tunneling currents. Finally, some applications for electron transfer reactions in models related to the bc1 complex of the respiratory chain will be presented.

PḊ = −

ℏ Im(Ψ*(t )∇⃗Ψ(t )) m

J ⃗ ·ds ⃗

(17)

Thus, in order to calculate the coupling, one needs to calculate the tunneling current J,⃗ integrate over the dividing surface, ∂ΩD, multiply by a fundamental constant ℏ, and add a minus sign. This general derivation can be repeated using a more general current density operator j ⃗ ̂ ( r ⃗) and multielectronic states |ΨD⟩ and |ΨA⟩, as in ref 5, instead of eq 9 that assumes explicitly the one-electron nature of the wave functions. In a more general

(8)

where VDA is the coupling (same as in eq 4) and ρ is the formal density of the Golden Rule states. (To avoid confusion, this density should not be identified with the thermally averaged 7713

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B

overlaps for such orbitals are close to unity, in practice no less than 0.1−0.5. The states therefore have the form

case, the transition current has the form of eq (2), everything else remains the same, including eq 17.

|D⟩ = |core(D), φtD⟩

3. ELECTRONIC FRANCK−CONDON FACTOR AND ONE TUNNELING ELECTRON APPROXIMATION According to eq 2, in the general multielectronic description, the transition flux (or current) J(⃗ r)⃗ is given by the matrix element of a one-electron operator j ⃗ ̂ ( r ⃗) between the donor and acceptor diabatic states |D⟩ and |A⟩. The two diabatic states |D⟩ and |A⟩ are nonorthogonal, as they belong to two different Hamiltonians (one in which the ground state is such that the transfer electron is on the donor and another in which the electron is on the acceptor). The calculation of such matrix elements is done most conveniently by using corresponding orbitals. Namely, when the donor and acceptor states are represented by the single determinant many-electron wave functions, a set of corresponding orbitals can be introduced by appropriate unitary transformation of the original HF canonical MOs of both states

|A⟩ = |core(A), φt A⟩

J ⃗ ( r ⃗) = (∏ sj) j≠t

OTE VDA = V DA FC(t )

FC(t ) =

∏ sj = ⟨core(A)|core(D)⟩ (25)

j≠t

OTE V DA =−

(18)

1 ⃗ J ( r ⃗) si i

ℏ (φ D∇φi A − φi A∇φi D) 2m i

(19)

(20)

∏ sj j≠i

ℏ2 2m

∫∂Ω

D

(φtD∇φt A − φt A∇φtD) d s ⃗

(26)

is a Koopmans-type, one tunneling electron (OTE) approximation for the coupling matrix element. Since the exchange effects are taken into account in the HF Hamiltonian, the OTE approximation can also describe the hole transfer.1b,6c,8 The form of the coupling matrix element in eq 24 is quite general; it suggests that VOTE DA can be calculated by a method not necessarily based on the tunneling currents. Indeed, the general form of the coupling matrix element eq 24 involving the electronic Franck−Condon factor FC was first described by Marshall Newton;9 he introduced the use of corresponding orbitals in the ET problem, and arrived at the concept of the electronic Franck−Condon factor in the analysis of the tunneling splitting based on the Fock operator. It is interesting to see how the same concept had resurfaced in the theory of tunneling currents; not surprising of course, because both theories express the coupling VDA as a matrix element between two nonorthogonal diabatic states of a one-electron operator (current density operator j ⃗ ̂ or Fock operator F̂). In the next section, we will expand on the equivalence of two approaches, and compare direct evaluation of the matrix element using Newton’s method and that based on tunneling currents. For proteins, the one-electron tunneling approximation holds particularly well;5b in this approximation, a single pair of the corresponding orbitals with the smallest overlap sσt gives the largest and dominant contribution to the tunneling current, whereas the rest (core orbitals) undergo induced polarization and participate as an electronic Franck−Condon factor.8 For example, in NADH dehydrogenase (complex I of the respiratory chain) for electron transfer between two FeS clusters N6a → N6b, the overlap of the tunneling orbitals is calculated to be sαt = 5.0 × 10−5. The next smallest overlap is 0.6 and all other overlaps are larger than 0.9, suggesting that core orbitals undergo only minor induced polarizations. The electronic Franck−Condon factors were typically in the range between 0.44 and 0.75 for all processes studied in complex I, where the donor and acceptor are multinuclear Fe/S clusters. Comparing the canonical and corresponding MOs, the atomic populations of the corresponding donor and acceptor

In this picture, each electron undergoes an independent transition from the initial to final state orbital, φDi → φAi , with corresponding flux given by eq 20. The total current is naturally the sum of currents between corresponding orbitals, one pair of orbitals per electron. Corresponding orbitals resemble normal modes of molecular vibrations. They are also closely related to Dyson orbitals (see, e.g., ref 7 and the discussion below). Notice that each one-electron term Ji ⃗ ( r ⃗) has a weight factor, which is the Franck−Condon overlap of the initial and final states, in which one pair of orbitals (and one electron) is excluded FC(i) = ⟨A|D⟩(i) =

(24)

where

so that the overlap matrix of |ΨD⟩ and |ΨA⟩ is diagonal: = δijsi.6 In the basis of corresponding molecular orbitals, the matrix element for current assumes a particularly simple form:

Ji ⃗ ( r ⃗) =

(23)

Hence, the coupling matrix element takes the following form:

⟨φAi |φDj ⟩

i

ℏ (φ D∇φt A − φt A∇φtD) 2m t

is the electronic Franck−Condon factor, and

|ΨA⟩ = |φ1A⋯φNA⟩

j

φAtun)

and a single pair of tunneling orbitals contributes mostly to total current. This allows introducing an effective one-electron picture of the transition:

|ΨD⟩ = |φ1D⋯φND⟩

J ⃗ ( r ⃗) = (∏ sj) ∑

(22)

(φDtun,

(21)

where sj = ⟨φAj |φDj ⟩. Naturally, in a tunneling transition, all electrons undergo a change of state, even if only a single charge is transferred; however, typically only a few pairs of orbitals with smallest overlaps si contribute significantly to the total current. This drastically simplifies the multielectronic picture of the transition. In particular, for long-distance electron tunneling, only one orbital undergoes significant “relocation” in space, some 10−20 Å, from being localized on the donor to that on the acceptor, and naturally the overlap for such a pair is very small, typically in the range 10−4−10−6; only far tails of such orbitals overlap. All other orbitals remain in the same region of space, pinned down to the same bond or group of atoms, and only very little shift or deform in response to a transferred charge, so that the 7714

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B

and acceptor orbitals takes place. This electronic time scale is the traversal time of the tunneling barrier discussed extensively in the past by Landauer and Buttiker,13 and most recently by Nitzan and co-workers in the context of electron transfer reactions.14 For our problem, this time scale can be estimated as

tunneling orbitals are significantly more localized on a single atom within the redox complex, suggesting that the electron tunneling occurs primarily between a pair of donor and acceptor atoms in the corresponding orbital picture. Concluding this section, we should point out that the tunneling orbitals are closely related to Dyson orbitals (see, e.g., ref 10). The relation is as follows φtD(1) = φi A(1) =

N ∏j ≠ t sj N ∏j ≠ i sj

τtun = Ntun

∫ ψND(1, ..., N )ψN*A−1(t)(2, ..., N ) d(2), ..., d(N ) ∫ ψNA(1, ..., N )ψN*D−1(t)(2, ..., N ) d(2), ..., d(N )

φAt

where and are the orbitals from which an electron leaves the donor complex and arrives to an acceptor complex, respectively, ΨN and ΨN−1 are the wave functions with N and (N − 1) electrons. This connection suggests a method of extending the description beyond the HF single determinant representation of diabatic donor and acceptor states adopted previously, for the multielectronic wave functions under the integrals of eq 27 should not necessarily be of HF type. The limitations of the HF approach are discussed next.

4. LIMITATIONS OF THE HARTREE−FOCK APPROXIMATION. ELECTRONIC REORGANIZATION ENERGY The form of the tunneling matrix element in eq 24 assumes a specific relation between the time scale of the tunneling electron and that of the core electrons;5b however, the exact criterion of applicability of this approximation and its possible extensions are not clear. In the absence of rigorous theory, an analogy can be a useful guide. It is interesting to notice that a similar form appears in proton coupled electron transfer (PCET) reactions,11 where the coupling has the form el VDA = V DA ⟨χ A |χ D ⟩

(28)

where is pure electronic coupling, and ⟨χ |χ ⟩ is the Franck−Condon factor for the initial and final vibrational states of the proton. Here the tunneling electron is coupled to a proton that changes its state (together with an electron) from the initial state χD to the final state χA. In PCET reactions, it is known that this is only one limit of the vibronic coupling which holds only for a weak electronic coupling. The applicability of eq 28 was carefully analyzed,11,12 and it was shown that it is related to a clear separation of proton tunneling time scale and electronic time scales. In general, one needs the proton time scale to be much shorter than that of electron transfer. Using this analogy, one can postulate that the applicability criterion of eq 24 should require the time scale of the tunneling electron to be much longer than that of the core electrons VelDA

τtunωcore ≫ 1

(30)

where we have assumed that the virtual orbitals making up the tunneling barrier are of the same energy as that of the bound core electrons, ΔEel, and Ntun is the number of tunneling steps (bridging atoms) that the electron is making to reach the acceptor, which is defined by the distance between the donor and the acceptor. The typical number of steps for a long distance tunneling would be 10 or more, Ntun > 10; thus, for a long distance tunneling, the criterion (29) appears to be satisfied. However, the issue is more complicated, as it also involves implicit use of a Born−Oppenheimer (BO)-type approximation1b,2 and its breakdown at the far tails of the tunneling electronic wave functions;15 here the role of the nuclear vibrations is played by the core electrons with their effective “vibrational” frequency ωcore ∼ ΔEel/ℏ, and the role of electronic wave function by the tunneling orbital that describes the tunneling electron. It is the tails of tunneling orbitals that define the coupling of redox centers and hence the rates of electron transfer at long distances. The Born−Oppenheimer approximation is equivalent to the Hartree−Fock selfconsistent picture of the electronic state, and its breakdown indicates the need of using post HF methods for the electron tunneling problem (i.e., extending electronic wave functions beyond the single determinant form, eq 18).16 To be specific, in the usual HF (or equivalently “BO”) calculations, the form of the wave function of, say, the donor state has the form of eq 22; i.e., the core electrons occupy orbitals that correspond to the tunneling electron being present at the donor site. However, the far tail of the tunneling orbital φDtun( r ⃗) corresponds more to a state of the core in which the tunneling electron is already withdrawn from the donor site, and therefore, the calculation of the tunneling tail of the donor/ acceptor orbitals should also involve a new state

(27)

φDt

ℏ ΔEel

A D

|D̃ ⟩ = |φt̃ D , core(A)⟩

(31)

(Or, in general, the intermediate between the D and A states of the core; see section 5.2.4.) The difference between states (22) and (31) is that in state (31) the tunneling electron has higher tunneling energy (and hence lower tunneling barrier) by the amount of the energy of reorganization of the core electrons from the donor to acceptor configuration. Since the overlap of the core orbitals is close to unity, one can assume that the difference should not be great; however, the issue has never been examined quantitatively and requires a separate detailed study. A quantitative measure of the effect is the energy of reorganization of the core electrons:

(29)

where τtun is the time scale of the tunneling electron, and ωcore is the frequency corresponding to core electrons; the latter is of the order of ΔEel/ℏ, and ΔEel ∼ 1 eV. Since in a HF selfconsistent approximation (eq 18) all electrons are treated on an equal basis, criterion (29) appears to be impossible to satisfy.5b In fact, more careful analysis shows that the time scale τtun is not the same as that of the bound motion of core electrons but is rather related to a nontrivial electronic time scale of reaching the far tail of the tunneling orbital, where the overlap of donor

λelD = E(φtD , core(A)) − E(φtD , core(D)) λelA = E(φt A , core(D)) − E(φt A , core(A))

(32)

When the core system does not change in the transition (as in the frozen-core approximation), the reorganization energy is 7715

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B

electrons and pseudopotentials for core electrons. The model system was fragmented into donor, bridge, and acceptor groups of different charge and spin multiplicity, and by this method, the tunneling electron was localized either on the donor or the acceptor sites (Figure S2, Supporting Information). The generated initial MOs were used in a subsequent BSZINDO calculation to find the donor and acceptor ground diabatic states,

close to zero. If the tunneling jump occurs quickly, so that τtunωcore ≪ 1, the core system does not have enough time to reorganize during the tunneling jump; however, the core gets excited upon the transition acquiring an additional energy λel (because after the transition, the core is still in the D state, while the “correct” state of the core is that of A). On the other hand, if transition occurs slowly, τtunωcore ≫ 1, the core adiabatically changes its state from D to A, always remaining in the ground state and acquiring no additional energy; in this case, all or part of the reorganization energy is transferred to the tunneling electron, thus reducing the tunneling barrier. We expect in this case significant modification of the tail of the tunneling orbital and much stronger coupling. The coupling of the tunneling electron to the remaining core electrons is similar to the vibronic coupling15c (see also the above-mentioned forthcoming paper on this subject by one of the authors), and it is clear that generally the problem amounts to dynamic correlation between electrons and requires methods beyond HF approximation. An interesting problem almost identical to that discussed in this paper involves a tunneling electron and image charge in solid state tunnel junctions.17 Here a tunneling electron is coupled to image charge, due to electron polarization of the electrodes, which plays the role of core electrons in this paper. An interesting theory was developed in ref 18, and experiment17 clearly indicated dynamic polarization effects (i.e., polarization that changes as the electron crosses the tunneling barrier), with a characteristic switch from BO to non-BO behavior at large tunneling distances. In these systems, the tunneling barriers are well controlled experimentally and can be up to 100 nm wide. In calculations described in the next section, we will mainly use the usual HF states but will also report some preliminary data on the extension of the HF method.

|ΨD⟩ = |φ1,Dσ , ..., φND, σ ⟩,

|ΨA⟩ = |φ1,Aσ , ..., φNA, σ ⟩

(33)

The obtained occupied canonical donor/acceptor MOs were then biorthogonalized2,8 to obtain the corresponding tunneling donor and acceptor orbitals; one pair of such orbitals shows the smallest overlap, while the other corresponding orbitals have overlaps close to unity: |Ψ′D⟩ = |ξ1,Dα , ..., ξl D, α , ξ1,Dβ , ..., ξtD, β⟩, |Ψ′A⟩ = |ξ1,Aα , ..., ξl A, α , ξ1,Aβ , ..., ξtA, β⟩

(34)

⟨ξi D, σ |ξjA, σ ⟩ = δijsiσ

(35)

sσi

The product of the overlap of the biorthogonalized states is the electronic Franck−Condon factor ∏li sαi ∏tj≠t sβj , where we assumed that the tunneling electron has a β-spin, and the tunneling orbital is the last one, having index “t” among the βspin corresponding orbitals. Tunneling matrix elements were calculated using interatomic tunneling currents. Since Gaussian ZINDO canonical MOs are represented in a delocalized basis set, obtained by orthogonalization23 of the original atomic orbitals, for proper treatment of interatomic tunneling currents, one needs to reverse orthogonalization and return to the original localized atomic orbitals. In localized but nonorthogonal atomic basis set, the interatomic currents have the following form:5

5. APPLICATIONS FOR ELECTRON TRANSFER IN A MODEL SYSTEM 5.1. Methods. 5.1.1. Model System Preparation. The model system was built based on the lower potential heme and higher potential heme (heme bL)/(heme bH) redox pair in the ubiquinol:cytochrome c oxidoreductase complex (known also as the bc1 complex or respiratory complex III) (PDB: 1BCC).19 The two heme molecules were extracted from the protein structure, along with their ligands, and the broken bonds were capped with the NHCH3 group for C-termini and the COCH3 group for N-termini; the heme propionate groups were kept neutral. To further simplify the system, the linker chain of repeating methylene groups was created manually and the model system was geometry-optimized using the PM6 semiempirical method implemented in the Gaussian 09 package,20 where only the heavy atoms of the heme molecules were kept fixed in space. The linkers of different lengths were generated by removing or adding one or two linker groups. Thereby, we have constructed five model systems M1−M5 with the same donor and acceptor groups and different linker lengths (Figure S1, Supporting Information). The models serve as a starting point of our upcoming systematic study of electron transfer in the bc1 complex. 5.1.2. Tunneling Current Calculation. The initial guess molecular orbitals (MOs) were generated by using the brokensymmetry (BS) B3LYP method21 in the Gaussian 09 package, using the basis set of the ZINDO method22 for the valence

Jab =

∏ siα ∏ siβ ∑ ∑ [(H̅ νμ − E0Sνμ)(θμν − θνμ) i



i≠t

ν∈a μ∈b

∑ (⟨νλ|μρ⟩ + SνμF ̅ )[(θμρBλν̅ β − θνρBλμ̅ β ) λρ

β β β β β β + (θλνBμρ ̅ − θλμBνρ ̅ )] + Sνμ ∑ [(Bνρ ̅ Bλμ ̅ − Bμρ ̅ Bλν ̅ ) λρ

× (∑ ⟨λρ|δγ ⟩θγδ)] (36)

γδ

Here ν and μ are the atomic orbitals of atoms a and b; θνμ = AνDμ, where Dμ and Aν are the expansion coefficients of the donor and acceptor tunneling orbitals, respectively; H̅ νμ and Sνμ are the reduced core Hamiltonian and the overlap matrix of the basis set; and E0 is a tunneling energy defined by E0 =

∑ DλFλρ̅ Dρ = ∑ AλFλρ̅ Aρ λ,ρ

(37)

λ,ρ

B̅ σλν

σ ∑i CσA iλ (1/si )

where F̅λρ is the reduced Fock matrix, and = CσD is the so-called reduced density matrix for the core orbitals, iν where CσA/D are expansion coefficients of donor/acceptor iλ molecular orbitals i with a given spin projection σ in atomic basis set λ, ν. The second equality in eq 37 corresponds to the resonance of the donor and acceptor potentials at the transition state of electron transfer reaction; in practice, the resonance is achieved by applying a static electric field mimicking the action 7716

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B

5.1.3. Complete Active Space Self Consistent Field (CASSCF) Calculation. To compare different schemes of electronic coupling calculation, on some models, we performed the CASSCF calculations which evaluated the avoided crossing splitting of the donor and acceptor states directly. To this end, the BS-ZINDO donor ground state was utilized; the delocalized canonical MOs were transformed into the natural orbitals (NOs) and they were used to generate the active space, which comprised six highest-occupied NOs. Three of them were localized on the donor iron atom, and all three were doubly occupied. The remaining three were localized on the acceptor iron atom, where only two were doubly occupied and the third was singly occupied by a α-spin valence electron. The search of the avoided crossing point was done by using Gaussian routines. This method was limited to systems with shortest links and hence largest splitting. 5.1.4. Direct Evaluation of the Coupling Using Corresponding Orbitals. In addition to tunneling currents and CASSCF, we applied the direct method of evaluating the coupling matrix element using corresponding orbitals.6c,24 In this method,

of the polar environment and solvation effects (Figure S3, Supporting Information). A detailed derivation of a somewhat lengthy formula (eq 36) is given in ref 5, and it involves a transition from continuously distributed in space flux to a coarse-grained picture of interatomic currents. The derivation and the final result become somewhat complicated due to nonorthogonality of the atomic basis functions on neighboring atoms, and the natural requirement of localization of the atomic orbitals to a given atom, so that individual atoms could be identified by their atomic basis functions. (The orthogonalization of the atomic basis sets does not help, as it results in the delocalization of the atomic basis functions, which one wants to avoid in the interatomic currents description.) For other details, one is referred to ref 5. The total atomic current through a given atom “a” is expressed as follows: Jatot ≡

1 2

∑ |Ja ,b | b

(38)

The total atomic current is proportional to the probability that the tunneling electron in passing through this atom; as such, it provides a convenient way (e.g., with color coding, Figure 1) of identifying atoms that constitute the tunneling pathways.

1 (E1 − E2) 2 ⎛ ⎞ 1 = ⎜HDA − SDA(HDD + HAA)⎟ /(1 − SDA 2) ⎝ ⎠ 2

VDA =

(39)

Here E1 and E2 are the energies of adiabatic states at the avoided-crossing point, SDA = ⟨D|A⟩ is the overlap between the donor and acceptor diabatic states, and HIJ are matrix elements of the Hamiltonian at the crossing point. In the Hartree−Fock approximation, using corresponding orbitals, the coupling takes the following form: VDA =

N ⎞⎤ det(S) ⎡ ⎛ 1 ̂ DA 1 ̂ DD ̂ AA )⎟⎥ ⎢∑ ⎜ f − + ( f f ⎟⎥ ii 2 ii 1 − det(S)2 ⎢⎣ i ⎝ si ii ⎠⎦ N



N

DA

∑ ∏ sj[fiî i

j≠i



AA 1 ̂ DD si(fii + fiî )] 2

N

=

∑ FC(i)ViiMO i

(40)

where fiĵ are matrix elements of the Fock operator, (S)ij = δijsi, and det(S) = ∏Ni si ≪ 1. The above expression contains the summation over N occupied corresponding MOs of tunneling matrix elements ViiMO multiplied by the corresponding electronic Franck−Condon factors FC(i). The one-electron approximation is derived from the above expression (cf. ref 9), where only one term with the largest FC is retained; all other terms are negligible, since the corresponding FC factor will be numerically very small. The convergence of the sum in eq 40 was numerically tested, as discussed in the next section. It has to be pointed out that, in addition to the selected methods reviewed above, there is significant literature on the coupling matrix element calculations; by no means do the reviewed methods cover the work of the past two decades on the subject, see, e.g., recent papers.25 5.2. Results and Discussion. 5.2.1. Tunneling Currents Method. Figures 1 and 2 and Table 1 summarize results for tunneling current calculations. The donor and acceptor diabatic

Figure 1. Total tunneling flux in model M1 across the dividing surface between the donor and acceptor redox sites. The red color intensity indicates the relative importance of the atoms in the tunneling process. The tunneling flux is expected to be relatively constant in the region between the two redox sites.

The interatomic currents are real, by their meaning, but they can be positive or negative, reflecting the direction of the current flow. In the above equation, we take the total value of the current, neglecting the sign. The difference in signs of currents ultimately reflects the quantum phase of the wave functions. The addition of currents of different signs represents the quantum interference effect, and as such reflects the true quantum nature of the flux. The calculation of the coupling matrix element is based on the flux theorem discussed in earlier sections, which for interatomic currents takes the form of eq 3. 7717

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B

Figure 2. Electronic coupling in models M1−M5 obtained by interatomic tunneling current theory TDA and by the direct calculation method, TMO DA , based on eq 40.

Table 1. Tunneling Calculation Results for the Five Model Systems Using the Interatomic Tunneling Current Method model M1 M2 M3 M4 M5

C-to-C distance (Å)

E-to-E distance (Å)

20.012 19.051 16.814 15.040 13.812

11.524 10.524 8.324 6.522 5.324

α-eFC

β-eFC

total eFC

⟨flux⟩ (cm−1)

⟨TDA⟩ (cm−1)

0.1877 0.5885 0.6525 0.8842 0.9064

0.8484 0.7022 0.8470 0.8707 0.8972

0.1592 0.4133 0.5527 0.7699 0.8132

1.427 2.569 1.44 × 101 5.07 × 101 1.26 × 102

2.272 × 10−1 1.0618 7.9630 3.901 × 101 1.025 × 102

TMO overlap 4.22 6.20 2.66 1.00 3.20

× × × × ×

−6

10 10−6 10−5 10−4 10−4

By applying an external electric field along the donor− acceptor axis, we bring the energy of the tunneling pair to resonance. Once at resonance, the tunneling flux across a separating surface was calculated using the interatomic tunneling currents. Results are shown in Figure 1. The degree of atomic contribution to the tunneling pathway is represented by the intensity of the red color. The tunneling current theory results are collected in Table 1. As expected, the tunneling flux and the corresponding tunneling matrix element drop exponentially as a function of the donor/acceptor distance, Figure 2. It is worth noticing that the electronic Franck− Condon factor in Table 1 decreases for larger systems. 5.2.2. Multielectronic States Avoided Crossing. Using the CASSCF method, the energy splitting between the donor and acceptor multielectronic adiabatic states was directly evaluated for two models, M4 and M5. Figure S4 (Supporting Information) shows the pattern of avoided crossing achieved by application of the external field. Table 3 summarizes results obtained with different methods, including CASSCF; the last column shows the calculated energy splitting (for additional details, see Table S1, Supporting Information); as expected, it decreases with the distance. The method is limited to systems where the energy splitting is sufficiently large; for this reason, we were unable to use it for systems with the longest links and small couplings. 5.2.3. Direct Method. Table 4 and Figures 2 and 4 summarize the results of the direct method evaluation of the tunneling matrix element for models M1−M5. As displayed in Figure 4, we have computed all partial contributions of terms in eq 40, starting from the tunneling individual VMO ii

ground states were calculated as discussed previously, and the biorthogonalization scheme was applied to find the corresponding orbitals. All corresponding core β-MOs show approximately unity overlap, except for the tunneling β-MO pair, which shows an overlap of order (10−4−10−6). The α-MOs exhibit close to unity overlap. The corresponding FC factors are shown in Table 1. The corresponding tunneling pair is a linear combination of all the canonical β-MOs, with significant contribution of many canonical β-MOs, Figure 3 and Table 2.

Figure 3. Composition of the donor and acceptor tunneling orbitals in terms of canonical MOs for model M1. For clarity, negligible overlaps with the low-lying canonical MOs are omitted from the above plots.

Table 2. Overlap Data between the Most Contributing Canonical β-MOs and the Donor/Acceptor Tunneling MOs (TMO) canonical HOMO TMO

−0

−2

−3

−4

−5

−7

−9

−10

−11

donor acceptor

−6.65 × 10−3 −0.929

−0.576 1.73 × 10−2

0.246 2.13 × 10−5

−0.716 3.26 × 10−5

−4.14 × 10−2 −0.279

−0.132 −3.10 × 10−2

0.147 −3.32 × 10−2

−0.116 4.17 × 10−2

−0.127 −5.32 × 10−2

7718

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B Table 3. Tunneling Calculation Results of Different Methods for All Five Model Systems ⟨TDA⟩ (cm−1)

model

−1 ⟨TMO DAtt⟩ (cm )

−1

2.272 × 10 1.0618 7.9630 3.901 × 101 1.025 × 102

M1 M2 M3 M4 M5

3.225 7.607 5.252 3.369 9.516

−1 converged ⟨TMO DA ⟩ (cm )

−1

× 10 × 10−1

3.521 8.240 5.590 3.567 9.791

× 101 × 101

(1/2)(EES − EGS) (cm−1)

−1

× 10 × 10−1 × 101 × 101

25.40 118.27

Table 4. Tunneling Calculation Results for the Five Model Systems Using the Direct Calculation Method model M1 M2 M3 M4 M5

̂ fDD tt (a.u.) −5.46 −5.48 −5.55 −5.64 −5.69

× × × × ×

̂ fAA tt (a.u.) 1

10 101 101 101 101

−5.37 −5.41 −5.49 −5.59 −5.67

× × × × ×

̂ fDA tt (a.u.) 1

10 101 101 101 101

−3.12 −3.29 −1.43 −5.44 −1.03

× × × × ×

−1 VMO tt (cm ) −4

10 10−4 10−3 10−3 10−2

2.025 1.841 9.505 4.376 × 101 1.170 × 102

−1 ⟨TMO DAtt⟩ (cm )

3.225 7.607 5.252 3.369 9.516

converged VMO (cm−1)

−1

× 10 × 10−1 × 101 × 101

−1 converged ⟨TMO DA ⟩ (cm )

2.211 1.994 1.011 × 101 4.633 × 101 1.204 × 102

3.521 8.240 5.590 3.567 9.791

× 10−1 × 10−1 × 101 × 101

core orbitals were formed as mixtures of donor and acceptor states in the from core(D̃ ) = (1 − X )core(D)0 + X core(A)0 core(Ã ) = (1 − X )core(A)0 + X core(D)0

(41)

With appropriate normalization, X represents the degree of mixing, and the “reaction coordinate” along which the core is changing from the initial (X = 0) to final (X = 1) state. The states used for calculation of the tunneling matrix element now have the form |D̃ ⟩ = |core(D̃ ), φtun ̃D ⟩

Figure 4. Partial sums over corresponding MOs in eq 40. The converged VMO is very close to that obtained for a single pair of tunneling orbitals represented by the first data point.

|Ã ⟩ = |core(Ã ), φtun ̃A ⟩

(42)

φ̃ D/A tun

where is the tunneling orbital reoptimized for a new core. Figure 5 show the results of calculation of the tunneling matrix

corresponding orbitals, down to the core ones. The data shown are for model M1 (the longest chain); the sum converges to a value of 2.21 cm−1, whereas the first term due to tunneling −1 orbitals VMO tt is 2.025 cm . This confirms our assumption that the main contribution comes from a single pair of tunneling orbitals. As shown in Table 4, which summarizes different contributions in eq 40, for all models, the converged ∑Ni VMO ii values are very close to those obtained with only one pair of tunneling orbitals VMO tt . The same trend holds for terms that MO include electronic FC factors, TMO DAi = FC(i)Vii . The distance dependence of the direct method evaluation is shown in Figure 2, where it is seen that the results are almost identical to those obtained with tunneling currents. 5.2.4. Extending the Hartree−Fock Picture. Following the discussion in section 4, here we explore an extension of the HF approach. To this end, we construct refined donor and acceptor states that should better represent the transition state of the system; namely, instead of the usual states represented by eq 22, we construct states in which the core orbitals are certain mixtures of the initial and final states, reflecting the gradual change of the core during the transition, and the tunneling orbital is readjusted to these new states of the core. As was discussed in section 4, such states, although still described by single determinants, phenomenologically should better represent the actual physical state of the electronic system at the transition state of the reaction. The new states by which the electronic coupling matrix element is now calculated were constructed as follows. The

Figure 5. Tunneling flux as a function of donor/acceptor (D/A) composition of the core. The solid line is the actual calculation; the broken line is the linear extrapolation to high values of D/A mixing. The optimization method becomes unstable at D/A composition X > 70/30, see text, as evidenced by the diverging solid line. At 50/50 mixing, the extrapolated tunneling flux is 6.812 cm−1, which corresponds to TDA = 1.085 cm−1. The HF value is T(HF) DA = 2.272 × 10−1 cm−1; see Table 1.

element using the new states. The tunneling matrix element, according to the flux theorem, is still a product of the flux, evaluated by using new φ̃ D/A tun , and the Franck−Condon factor. As seen in Figure 5, as we adjust the core, the tunneling orbitals modify, and initially the resulting flux almost linearly increases, reflecting the adjustment of the tails of the tunneling orbitals to the new state of the core; however, at some point 7719

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B (about X = 0.6), the optimization procedure that we use becomes unstable (it was no longer possible to converge the optimization at a given state of the core) and the resulting flux begins to markedly deviate from the initial behavior. At ratio 50/50, X = 0.5, the optimization procedure could not be used at all. (Our current procedure of reoptimization of tunneling orbitals with the fixed core is obviously not optimal, and the work is underway to improve the procedure.) To estimate the matrix element, we therefore used extrapolation of the initial trend, which allows one at least semiquantitatively to estimate the magnitude of the effect. As seen, the extrapolated value of the tunneling matrix element at 50/50 ratio, which corresponds to the transition state, is a factor of 3−5 different from the HF result; see the Figure 5 caption. Although a more detailed study of the effect is obviously needed, the preliminary results indicate a significant effect that can give an order of magnitude difference in the rate (which is the square of the coupling).

role, indicating explicit limitations of the HF approach. Although the physics of the limitation factors is clear, the computational methods are nontrivial and yet to be fully developed. Clearly, much work on the nature of donor/ acceptor interactions, pioneered by Marshall Newton and John Miller some three decades ago, is still needed.



ASSOCIATED CONTENT

S Supporting Information *

Results of CASSCF calculations for models M4 and M5, the reorganization energies calculated for the five models, the shortest and longest M1 and M5 models, atomic populations of donor/acceptor tunneling orbitals for M1 model, resonance energies of the donor and acceptor tunneling orbitals obtained by applying an external electric field, and the avoided crossing between the ground and first excited multi-electronic adiabatic states for model M5. This material is available free of charge via the Internet at http://pubs.acs.org.



6. SUMMARY AND CONCLUSIONS Table 3 summarizes all the results obtained for models M1− M5. As seen, for large and moderately small couplings, where direct methods can be applied, the tunneling currents reproduce results of other methods rather accurately. In other cases, when the coupling is extremely small, the direct methods either become impossible to apply or become unreliable, as it becomes difficult to assess their accuracy. The tunneling currents method, on the other hand, as it deals with the tails of wave functions directly and has an internal check of accuracy (the invariance of the flux calculated for different dividing surfaces), has the best chance to evaluate the tunneling coupling most accurately. Some examples of extremely small couplings evaluations (in the range of 10−2 cm−1 and below) were demonstrated in our previous publications, e.g., refs 1b and 3. In this paper, we have provided an elementary proof of the key result on which the tunneling current method is based, namely, the flux theorem, by using the framework of the Golden Rule. This approach is more natural, compared with what we had before, as it deals with the rate theory from the start, and shows explicitly how the tunneling currents appear naturally in the expression of the coupling matrix element. This is an important step in the electronic coupling theory, in particular as it has a natural extension to multielectronic formulation of the problem. The method however is limited by the use of the Hartree− Fock approximation. We have discussed the HF limitations in some detail in sections 4 and 5. The key issue is the correlation between the tunneling electron and the core electrons. An indirect quantitative measure of the effect is the electronic reorganization energy, eq 32, which amounts to a few eV; see Table S2 (Supporting Information). In order to directly probe the effect on the coupling, here we presented our first results with phenomenologically adjusted states that should better represent the transition state of the system. The preliminary results indicate that non-Hartree−Fock effects can amount to an order of magnitude change of the rate constant. What is important is that this conclusion should not depend on the method of calculation; indeed, here we showed that, in all cases where direct comparison was possible, the flux theory produces the same results as other methods. The difference is due to the nature of the diabatic donor and acceptor states used in the evaluation of the coupling. It appears from the present study that non-Hartree−Fock correlation effects play a significant

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS It is a great pleasure to contribute to this issue of The Journal of Physical Chemistry B, and to join other colleagues honoring two outstanding researchers in the field of electron transferDrs. Marshall D. Newton and John R. Miller. One of the authors (A.A.S.) wishes to express particular gratitude to Marshall Newton for many stimulating scientific discussions and friendship over the past two decades. Supported by NIH grant GM54052.



REFERENCES

(1) (a) Stuchebrukhov, A. A. Toward ab initio theory of longdistance electron tunneling in proteins: Tunneling currents approach. Adv. Chem. Phys. 2001, 118, 1−44. (b) Stuchebrukhov, A. A. Longdistance electron tunneling in proteins. Theor. Chem. Acc. 2003, 110, 291−306. (2) (a) Newton, M. D. Electronic Coupling in Electron Transfer and the Influence of Nuclear Modes: Theoretical and Computational Probes. Theor. Chem. Acc. 2003, 110, 307. (b) Newton, M. D.; Sutin, N. Electron transfer reactions in condensed phases. Annu. Rev. Phys. Chem. 1984, 35, 437−80. (3) Hayashi, T.; Stuchebrukhov, A. A. Electron tunneling in respiratory complex I. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 19157−19162. (4) Formally, as in the previous derivation eq 5, we require orthogonality of the donor and acceptor states here; however, this is not necessary for the final expression of the coupling. The argument can be repeated with the assumption of nonorthogonality as well, and the final expression is the same. (5) (a) Stuchebrukhov, A. A. Tunneling currents in long-distance electron transfer reactions. IV. Many-electron formulation. Nonorthogonal atomic basis sets and Mulliken population analysis. J. Chem. Phys. 1998, 108, 8510−8520. (b) Stuchebrukhov, A. A. Tunneling currents in long-distance electron transfer reactions. V. Effective one electron approximation. J. Chem. Phys. 2003, 118, 7898−7906. (6) (a) Cory, M. G.; Zerner, M. Projected unrestricted Hartree-Fock calculations within the intermediate neglect of differential overlap model. J. Phys. Chem. A 1999, 103, 7287−7293. (b) Amos, A. T.; Hall, G. G. Single Determinant Wave Functions. Proc. R. Soc. London, Ser. A 1961, 263, 483−493. (c) Newton, M. D. Quantum Chemical Probes 7720

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721

Article

The Journal of Physical Chemistry B of Electron-Transfer Kinetics - the Nature of Donor-Acceptor Interactions. Chem. Rev. 1991, 91, 767−792. (d) Goddard, W. A.; Ladner, R. C. Generalized Orbital Description of Reactions Small Molecules. J. Am. Chem. Soc. 1971, 93, 6750−6756. (e) Zhang, L. Y.; Friesner, R. A.; Murphy, R. B. Ab Initio Quantum Chemical Calculation of Electron Transfer Matrix Element for Large Molecules; Schrodinger, Inc.: Portland, OR, 1997. (7) Ortiz, J. V. Brueckner orbitals, Dyson orbitals, and correlation potentials. Int. J. Quantum Chem. 2004, 100, 1131−1135. (8) Newton, M. D.; Ohta, K.; Zhong, E. Analysis of Superexchange Coupling in Metallocene-Metallocenium Redox Pairs. J. Phys. Chem. 1991, 95, 2317−2326. (9) Newton, M. D. Electronic-Structure Analysis of Electron-Transfer Matrix-Elements for Transition-Metal Redox Pairs. J. Phys. Chem. 1988, 92, 3049−3056. (10) (a) Oana, C. M.; Krylov, A. I. Dyson orbitals for ionization from the ground and electronically excited states within equation-of-motion coupled-cluster formalism: Theory, implementation, and examples. J. Chem. Phys. 2007, 127;. (b) Ortiz, J. V. Ionization energies and Dyson orbitals of 1,2-dithiin. J. Phys. Chem. A 2002, 106, 5924−5927. (11) Hammes-Schiffer, S.; Stuchebrukhov, A. A. Theory of Coupled Electron and Proton Transfer Reactions. Chem. Rev. 2010, 110, 6939− 6960. (12) Georgievskii, Y.; Stuchebrukhov, A. A. Concerted electron and proton transfer: Transition from nonadiabatic to adiabatic proton tunneling. J. Chem. Phys. 2000, 113, 10438−10450. (13) (a) Landauer, R.; Martin, T. Barrier interaction time in tunneling. Rev. Mod. Phys. 1994, 66, 217. (b) Buttiker, M.; Landauer, R. Traversal time for tunneling. Phys. Rev. Lett. 1982, 49, 1739. (14) Nitzan, A.; Jortner, J.; Wilkie, J.; Burin, A. L.; Ratner, M. A. Tunneling Time for Electron Transfer Reactions. J. Phys. Chem. B 2000, 104, 5661−5665. (15) (a) Beratan, D. N.; Hopfield, J. J. Failure of the BornOppenheimer and Franck-Condon approximations for long distance electron tranfer rate calculations. J. Chem. Phys. 1984, 81, 5753−59. (b) Onuchic, J. N.; Beratan, D. N.; Hopfield, J. J. Some aspects of electron-transfer reaction dynamics. J. Phys. Chem. 1986, 90, 3707− 3721. (c) Medvedev, E. S.; Stuchebrukhov, A. A. Breakdown of the Born-Oppenheimer-Condon-Marcus approximation in long distance electron transfer. Chem. Phys. 2004, 296, 181−192. (16) An extended discussion of the failure of BO approximation in a framework suitable for the electronic problem generalization is presented in a forthcoming paper in this journal by one of the authors (A.A.S.). (17) Gueret, P.; Marclay, E.; Meier, H. Investigation of Possible Dynamic Polarization Effects on the Transmission Probability of NGaas/Alxga1-Xas/N-Gaas Tunnel Barriers. Solid State Commun. 1988, 68, 977−979. (18) Persson, B. N. J.; Baratoff, A. Self-Consistent Dynamic Image Potential in Tunneling. Phys. Rev. B 1988, 38, 9616−9627. (19) Zhang, Z.; Huang, L.; Shulmeister, V. M.; Chi, Y. I.; Kim, K. K.; Hung, L. W.; Crofts, A. R.; Berry, E. A.; Kim, S. H. Electron transfer by domain movement in cytochrome bc1. Nature 1998, 392, 677−84. (20) 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.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, N. J.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09; Gaussian, Inc.: Wallingford, CT, 2009.

(21) (a) Becke, A. D. Density-functional thermochemistry. III. The role of exact exchange. J. Chem. Phys. 1993, 98, 5648−5652. (b) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J. Ab Initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields. J. Phys. Chem. 1994, 98, 11623−11627. (22) Thompson, M. A.; Zerner, M. C. A theoretical examination of the electronic structure and spectroscopy of the photosynthetic reaction center from Rhodopseudomonas viridis. J. Am. Chem. Soc. 1991, 113, 8210−8215. (23) Löwdin, P.-O. A Quantum Mechanical Calculation of the Cohesive Energy, the Interionic Distance, and the Elastic Constants of Some Ionic Crystals. Ark. Mat. Astr. Fys. A 1947, 35, 9. (24) (a) Curtiss, L. A.; Naleway, C. A.; Miller, J. R. Theoretical-Study of Long-Distance Electronic Coupling in H2c(Ch2)N-2ch2 Chains, N = 3−16. J. Phys. Chem. 1993, 97, 4050−4058. (b) Curtiss, L. A.; Naleway, C. A.; Miller, J. R. Superexchange Pathway Calculation of Electronic Coupling through Cyclohexane Spacers. J. Phys. Chem. 1995, 99, 1182−1193. (25) (a) Beratan, D. N.; Onuchic, J. N. Redox redux. Phys. Chem. Chem. Phys. 2012, 14, 13728−13728. (b) Skourtis, S. S.; Waldeck, D. H.; Beratan, D. N. Fluctuations in Biological and Bioinspired ElectronTransfer Reactions. Annu. Rev. Phys. Chem. 2010, 61, 461−485. (c) Nishioka, H.; Ando, K. Pathway analysis of super-exchange electronic couplings in electron transfer reactions using a multiconfiguration self-consistent field method. Phys. Chem. Chem. Phys. 2011, 13, 7043−7059. (d) Nishioka, H.; Ando, K. Electronic coupling calculation and pathway analysis of electron transfer reaction using ab initio fragment-based method. I. FMO-LCMO approach. J. Chem. Phys. 2011, 134;. (e) Sumi, H.; Kakitani, T. Unified Theory on Rates for Electron Transfer Mediated (vol 105, pg 9603, 2001). J. Phys. Chem. B 2009, 113, 12852−12852.

7721

DOI: 10.1021/jp512699a J. Phys. Chem. B 2015, 119, 7712−7721