The treatment of ET reactions in the usual Marcus picture assumes a canonical equilibrium distribution of all possible microstates of donor and acceptor, as implied by the brackets ⟨···⟩. This was also a crucial assumption in the derivation of the SE, FR, and hopping models in section 5.1 . The assertion of equilibrium thermodynamics can be expected to give a good approximation if the ET rate is slower than the slowest system vibrational mode coupling to ET. (157) Molecular dynamics simulations have shown that for small proteins the slowest modes coupling to ET are due to the electroelastic fluctuations of the protein/water interface occurring on the 1–10 ns time scale. (36) Even longer time scales may be possible, e.g., when large protein conformational changes couple to ET (see section 5.3 ). Hence, one can expect that for slow to modestly fast biological ET reactions occurring on the millisecond to microsecond time scale, equilibrium statistical mechanics and Marcus theory give an adequate description. Indeed, the successful interpretation of many experimental measurements of millisecond to microsecond biological ET reactions are an undeniable proof for the strength of this theory.

However, some biological ET reactions are ultrafast, occurring on the nanosecond to picosecond (ns–ps) time scale, e.g., heme a to a 3 electron tunneling in cytochrome c oxidase (ns), hole transport in DNA, DNA photolyase, and cryptochrome (ns–ps), and primary charge separation in photosystem II (ps). In these systems the ET is comparable to or faster than the electroelastic fluctuations at the biomolecule/water interface. Hence, a significant fraction, if not all, of these fluctuations are frozen on the time scale of the ET, invalidating the use of canonical equilibrium statistical mechanics. Consequently, the time averages measured in experiment are no longer equal to the canonical equilibrium averages, which means that the system is nonergodic.

Application of eq 85 to the picosecond primary charge separation reaction in bacterial photosynthesis has given a spectacular reduction in reorganization free energy from unphysically high values of 2.36 eV obtained from integration over all frequencies ( eq 62 ) to about 0.36 eV (233, 234) in good agreement with the experimental estimate of 0.22 eV; (24) see Table 1 . A similar analysis for hemetotunneling in cytochromeoxidase is discussed in section 6.2 . In general, the result of dynamical freezing of the slow modes is that reorganization and activation free energy as well as the magnitude of the free energy difference decrease compared to canonical equilibrium conditions, as shown schematically in Figure 12 B.

Figure 12. (A) Self-consistent correction scheme for reorganization free energy for fast (nonergodic) ET reactions. The lower integration limit (2π k n –1 ) for the spectral density function J ( eq 29 ) is iterated until it is equal to the ET rate k . The result is a nonergodic reorganization free energy, λ( k ), that is smaller than the equilibrium reorganization free energy including all frequencies, λ. In (B) parabolic free energy profiles are drawn for the equilibrium case (blue) and after the nonergodicity corrections have been applied (red). The nonergodicity correction for the free energy difference between the two states, i.e., the vertical difference between the minima, is assumed to result in the same scaling factor as for reorganization free energy. The Stokes shift is indicated by twice the reorganization free energy. See section 5.2.1 for further details.

If one assumes that the ET rate constant is known from experiment, say, one can replace in the formula for the canonical reorganization free energy in the frequency domain, eq 62 , the lower integration limit “0” by “”. This leads to the definition of a noncanonical (or nonergodic) reorganization free energywhich includes only the contributions from modes that are at least as fast as the reaction rate. However, these fast contributions are averaged over all of configuration space including the part described by the slow modes, because the spectral density functionis formally obtained from infinitely long trajectories. A similar correction has been suggested for the rate dependent average energy gap Δ (36, 234) and for the free energy difference Δ) by virtue of Δ) = Δ) – λ() . In case the ET rate constant is not known experimentally, one could solve the equations for λ() and Δ) together with the one for the ET rate (, e.g., eq 25 ) self-consistently until; see scheme in Figure 12 A.

Matyushov and co-workers suggested a simple correction for nonergodic effects based on the theory of dynamically restricted canonical ensembles. (262) The formalism is similar to equilibrium statistical mechanics with the crucial difference that integration is not done over full phase space, but over the part of phase space that is accessible on the time scale of the ET event. As a consequence, reorganization free energy, average energy gap, free energy difference, and activation free energy become a function of the ET rate.

A note of caution may be appropriate at this point. In certain systems the nonergodicity correction may lead to small reorganization free energies on the order of the electronic coupling matrix element. Then one needs to verify if the assumption of localized charge carriers and the use of the nonadiabatic or adiabatic rate formalisms ( eqs 24 26 ) are still appropriate, or if the carrier is delocalized and explicit charge propagation schemes should be used, as discussed in section 5.2.3

Matyushov and co-workers have stressed that also in this formalism a nonergodicity correction is necessary in the form of a self-consistent ET rate-adjusted diffusion coefficient) and free energy surface,(Δ) →,Δ) . Indeed, the authors have shown that self-consistent solutions to eqs 87 and 88 using data from MD simulation give good agreement with the multiexponential decay kinetics reported for a number of mutants of (112) The dynamical quenching of the slow protein response was suggested to be an important mechanism enabling ultrafast biological ET while reducing the loss of free energy to heat. (36)

The elimination of slow modes according to the nonergodicity correction discussed in section 5.2.1 gives self-consistent rate-adjusted ET parameters. When used in the framework of Marcus theory, the kinetics described by these parameters is still monoexponential. More complicated multiexponential decay kinetics, as measured, e.g., for photosynthetic reaction center proteins, requires a treatment beyond the standard Marcus formalism. (112) Multiexponential decay implies that motions on at least two different time scales couple to ET. In their theory developed some time ago, Sumi and Marcus describe the ET by two coordinates, one coordinatefor a fast polarization response and another coordinatefor a slow polarization response compared to the ET time scale. (119) In case of the picosecond photosynthetic charge separation reactionwould correspond to the subpicosecond response of the intramolecular vibrations of the cofactors, andwould correspond to all other relaxation processes of the protein and solvent. The total energy gap is assumed to be a linear function ofandresulting in a fast energy gap component, Δ, and a slow component, Δ. The fast component Δis assumed to be in thermal equilibrium for each value of Δ, whereas the slow component Δis assumed to diffuse on a free energy surface(Δ) with diffusion coefficient. The rate for ET along the Δdirection is(Δ) . In this picture the initial state populationis a function of Δand time only,(Δ), with a time dependence given by the diffusion-reaction (Fokker–Planck) equation (119) and the ET rate is obtained from the initial state population decay

5.2.3 Nonadiabatic Molecular Dynamics

k) ≫ H ab . An interesting situation occurs when this is no longer true. In H ab . If λ/H ab ≤ 2, the activation barrier for ET vanishes according to A = 0), and the transition state becomes a minimum. In this regime the excess electron is no longer localized on one cofactor or molecular fragment, but spontaneously delocalizes; the ET rate becomes ill-defined. This scenario may occur for two or more closely spaced cofactors or DNA bases with strong electronic interactions and surrounded by a low dielectric and/or slowly responding (viscous) medium. In all previous sections it was assumed that initial and final ET states are characterized by potential wells that are sufficiently deep so that localized charge carriers can form in initial and final states. This is the case when reorganization free energy, after correction for possible nonergodic effects according to eq 85 , is still significantly larger than electronic coupling, λ() ≫. An interesting situation occurs when this is no longer true. In Figure 13 I show the adiabatic free energy curves for different values of the ratio λ/. If λ/≤ 2, the activation barrier for ET vanishes according to eqs 18 and 19 (at Δ= 0), and the transition state becomes a minimum. In this regime the excess electron is no longer localized on one cofactor or molecular fragment, but spontaneously delocalizes; the ET rate becomes ill-defined. This scenario may occur for two or more closely spaced cofactors or DNA bases with strong electronic interactions and surrounded by a low dielectric and/or slowly responding (viscous) medium.

Figure 13 Figure 13. ET free energy curves for large ratios H ab /λ. The diabatic free energy curves eqs 11 and 12 are drawn in blue for λ = 0.4 eV and ΔA = 0 eV. Adiabatic ground and excited state free energy curves are obtained according to eq 10, and drawn in red for increasing values of electronic coupling H ab at constant λ = 0.4 eV. At H ab = λ/20 the ET can be classified as nonadiabatic, and at λ/5 it can be classified as adiabatic. For H ab ≥ λ/2 the free energy barrier disappears and the transition state at ΔE = 0 becomes a minimum.

3 O 4 , tetrathiafulvalene The phenomenon of excess charge delocalization is well-known in the chemistry literature on mixed-valence compounds, where the degree of delocalization is classified according to Robin and Day. (263) In class I compounds the excess charge is highly localized, in class II compounds there is some localization of distinct valences but with low activation energy for interconversion, and in class III compounds the excess charge is completely delocalized. Textbook examples for class I, II, and III compounds are Pb, tetrathiafulvalene (264) (also correctly predicted by CDFT (130, 205) ) and the Creutz–Taube ion. (265)

Considering a chain of closely spaced ionizable sites, the situation also has some similarity with the band transport problem in solids. However, in contrast to solids, the thermal fluctuations in proteins are so large that the mean free path of the electron is smaller than the typical site spacings, making band-theoretical approaches inapplicable. In other words, in this regime the electronic interaction between the sites can no longer be treated as a perturbation as is done in Marcus theory, while the strong nuclear fluctuations in proteins cannot be treated as a perturbation to the electronic problem as is done in band theory. Nuclear and electronic motion can no longer be safely separated and the coupled nuclear–electronic problem should be solved instead.

Potential solutions to this problem are direct propagation schemes of the coupled electron–nuclear motions. In particular, mixed quantum–classical (MQC) nonadiabatic molecular dynamics (NAMD) methods have been developed for (bio)molecular systems including the very early work by Warshel and co-workers, (266, 267) with applications to charge separation in photosynthetic reaction center proteins, (268) mean field (MF) Ehrenfest, and fewest switches surface hopping (SH). (269-272) These approaches are also referred to as semiclassical trajectory methods because nuclear motion is treated using classical mechanics and electronic motion is treated quantum mechanically. MF and SH have served as the methods of choice for the last 20 years, primarily in the context of ab initio and DFT electronic structure calculations of photoexcitation processes. (272) Also electron transfer reactions in the gas phase have been investigated by combining the MF approach with real-time TDDFT. (273) (In the latter study the competition between superexchange and hopping was investigated for a small donor–bridge–acceptor system.) Recently, the MF and SH methods have been adopted for simulation of thermal charge transport in larger, biological systems such as DNA and proteins. (37, 38, 250, 274, 275) I note in passing that charge transfer in organic semiconducting materials is in a similar regime, (199, 276, 277) and first NAMD simulations for simple one-dimensional model Hamiltonians (278-281) and approximate NAMD simulations for a slab of organic molecules (282) have recently been reported.

t) is expanded in a set of adiabatic electronic states ϕ l (ϕ 0 and ϕ 1 for a simple two-state donor–acceptor problem) (89) c l (t) the time dependent expansion coefficients and R the position vector of all nuclear coordinates. Insertion of (90) H kl are the Hamiltonian matrix elements, d kl are the nonadiabatic coupling vectors between states k and l, and Ṙ is the velocity vector for all nuclei. The expansion coefficients are propagated in time according to E k = H kk (e.g., E 0 , E 1 given by c k |2; see E k at any time and transitions between surfaces E l ← E k occur stochastically according to a transition probability derived by Tully; MF and SH are MQC schemes where the electronic wave function Ψ() is expanded in a set of adiabatic electronic states ϕ(ϕand ϕfor a simple two-state donor–acceptor problem)with) the time dependent expansion coefficients andthe position vector of all nuclear coordinates. Insertion of eq 89 in the time-dependent electronic Schrödinger equation giveswhereare the Hamiltonian matrix elements,are the nonadiabatic coupling vectors between statesand, andis the velocity vector for all nuclei. The expansion coefficients are propagated in time according to eq 90 for a time dependent electronic Hamiltonian as determined by the classical nuclear motion. The crucial difference between MF and SH is the way how the classical nuclear dynamics is treated. In MF the nuclei move on a mean field electronic surface composed of the adiabatic potential energy surfaces(e.g.,given by eq 3 for a two-state system) with weights proportional to the electron coefficients |; see Figure 14 A. In the SH method the nuclei are propagated on a single adiabatic surfaceat any time and transitions between surfacesoccur stochastically according to a transition probability derived by Tully; (269) see Figure 14 B. In both cases initial state population decay occurs when the nuclear dynamics generates nuclear configurations where the initial electronic state becomes (quasi)-degenerate with an excited electronic state. In practice, one needs to run a large number of Ehrenfest of SH trajectories with different initial conditions. Then effective ET rates can be obtained from the ensemble averaged time-dependent decay of the initial charge population and charge mobilities calculated from the mean square displacement of the center of excess charge versus time.

Figure 14 Figure 14. Two popular quantum–classical nonadiabatic molecular dynamics (NAMD) methods: (a) mean-field (MF) Ehrenfest MD and (b) surface hopping (SH) MD. For each case, two adiabatic potential energy surfaces are shown, E 0 and E 1 . The trajectories of two independent simulations are indicated in blue. After passing the high-coupling region at the avoided crossing, the trajectories evolve on a potential that is an average of E 0 and E 1 in the MF simulations. By contrast, the trajectories evolve on either E 0 or E 1 in the SH simulations. Reprinted with permission from ref 38. Copyright 2013 The Royal Society.

The strength of NAMD approaches is that the ET mechanism is, in principle, a result of the calculation, in contrast to Marcus or band theory, which are based on a preconceived picture of the ET process. However, the computationally amenable MF and SH approaches still have a number of well-known shortcomings (see below) and a major drawback is their very high computational demand, in particular when the electronic structure problem is solved at the ab initio (283) or DFT level. (273, 284) Elstner and co-workers have addressed the problem of efficiency by developing MF and SH implementations based on self-consistent charge density functional tight binding (SCC-DFTB) calculations of the electronic Hamiltonian. This allows for an efficient propagation of the coupled electron–nuclear motions for relevant biological systems on the picosecond–nanosecond time scale. I refer here to two recent reviews by the authors for a detailed presentation of their implementation. (37, 38) Applications of the method to DNA and DNA photolyase are discussed in sections 6.4 and 6.5

The shortcomings of the MF approach are well documented (see ref 285 and references therein): (1) while a reasonably good approximation to the short time dynamics of the system is obtained, after reaching an avoided crossing the system remains in a mixed state indefinitely instead of returning to one of the adiabatic states. A consequence of this is unphysically large delocalization of the excess electron or hole. (2) The method does not satisfy a detailed balance between quantum and classical subsystems, resulting in the quantum system acquiring too much energy during the course of the reaction. (3) A quantum–classical method, the MF method neglects nuclear quantum effects, which can become important at low temperatures. Issue 1 has been addressed by Jasper and Truhlar in their decay-of-mixing approaches, which impose demixing of the electronic states, i.e., a decoherence correction, to model the dissipative effect of the environment on the electronic dynamics. (286-288)

H ab ) and/or strong nonadiabatic coupling (d kl ). Recently, nuclear quantum effects were incorporated in the SH method The SH method and its latest extensions and modifications also improve on these issues. The stochastic hops mimic the effect of branching of the nuclear wave packets at the crossing region and the nuclei are propagated at any time on a pure adiabatic state. In recent analyses of the equilibrium limits of the SH method for a two-level (289) and a three-level model system, (290) it was found that SH does not in general yield a Boltzmann equilibrium population of the electronic states but that in practice the observed deviations are small, especially in the limits of small electronic coupling () and/or strong nonadiabatic coupling (). Recently, nuclear quantum effects were incorporated in the SH method (291) by combining it with ring polymer molecular dynamics (39, 292) (see also section 5.4.2 ). The tunneling contribution to the rate for a one-dimensional two state model of a nonadiabatic reaction could be well reproduced when compared with exact quantum mechanical calculations.

A longstanding problem of standard SH that is particularly relevant for ET simulations is the missing decoherence of the electronic wave function. After passing the crossing region, the off-diagonal electronic density matrix element does not decay sufficiently, leading to an overly coherent electronic wave function with probability density on two or more surfaces. This shortcoming is not new and has been addressed by many researchers. (286, 293-297) One of the simplest approaches is to collapse the electronic wave function to a single adiabatic state after leaving the crossing region. (294) More recently, the issue was readdressed by Subotnik and co-workers in a series of papers. Using a simple spin-boson model, the authors found that standard SH rates without decoherence correction scale incorrectly with electronic coupling. (298) They suggested an augmented SH algorithm where decoherence is introduced by stochastically collapsing the electronic wave function according to a rate that depends on positions and momenta. (296, 297) The new algorithm, denoted A-FSSH, was shown to recover the correct scaling of ET rates with electronic coupling. (297)