New Hamiltonian model for long-range electronic superexchange in

J. Phys. Chem. 1993,97, 5581-5593

5581

New Hamiltonian Model for Long-Range Electronic Superexchange in Complex Molecular Structures James M. Gruschus and Atsuo Kuki’ Department of Chemistry, Baker Laboratory, Cornell University. Zthaca, New York 14853- I301 Received: January 26, I993

The electronic superexchange interactions, which enable long-range electron transfer in complex molecular structures, such as proteins, are beyond the capacity of standard electronic structure methods due to their size, inhomogeneity, aperiodicity, and sensitivity to the many weak interactions between the nominally insulating atoms of the intervening medium. The inhomogeneous aperiodic lattice theory, presented here, is a novel strategy implemented as a quantum molecular model, designed to enable the calculation of electronic transfer matrix elements in large macromolecular systems by redesigning the electronic structure problem into a twotiered approach. The procedure is (i) assembly of the diagonal and off-diagonal elements of the Hamiltonian matrix by comparison, respectively, with experimental ionization potentials for individual amino acids and with triple-cab initio studies of resonance integrals and then (ii) nonperturbative computation of the macromolecular electronic coupling from this inhomogeneous aperiodic lattice (IAL) Hamiltonian. The IAL method includes all the occupied orbitals of the entire protein in a full-matrix calculation with no a priori assumption about pathways or spatial subregions important in spanning the distance between the redox sites. The nonperturbative approach and overall strategy of subunit-level calibration, distinct from standard semiempirical atom-level calibration, are first presented. A specific algorithm for Hamiltonian construction and tests against a series of medium sized molecules then follow. This nonperturbative matrix inversion strategy is computationally efficient and can treat a system the size of cytochrome c in less than a minute. Most important to the mechanistic study of redox proteins, the absolute results in cm-l for the computed charge resonance energies of three ruthenium modified cytochrome c derivatives compare well with experiment.

with carefully calibratedproperties. The properties of the subunits themselves are not computed but rather are drawn from small Recent experimentson the intramolecular aspects of electronmolecule experiment (on noninteracting single moleculcs) and transfer kinetics in respiratoryand photosynthetic proteins, which ab initio studies (on interacting small molecule pairs). Thus the are central processes in biological energy transduction, have key idea is a strategy for macromolecular model Hamiltonian evolved to an advanced level due particularly to the availability construction, which is in a sense intermediatebetween a continuum of X-ray crystal structures and the concomitant development of model, based on exponential and isotropic falloff through some semisyntheticmodification and site-directed mutagenesis capaeffective barrier, and a standard molecular electronic structure bilities. I-3 These experimentshave made increasingly clear that perspective, in which the aim may be to compute SCF molecular electron transfer can proceed surprisingly rapidly as a single orbitals for all electronsor for all valence electrons. Instead, this elementary step over long range: and this phenomenon has is a twetieredstrategy,which does not attempt to calculate small stimulated theoretical investigations into fundamental aspectsof molecule and macromolecular electronic structure simultaneously, the quantum electronics of proteins. Perturbation theoretical but rather (i) installs realistic molecular component properties methods are well-established in the treatment of the primary from reliable sources and then (ii) uses nonperturbativenumerical process in bacterial photosynthesis,5.6but new and distinct features Green’s function evaluation to calculate the macromolecular arise at longer range where many intermediary amino acids must coupling. In this first paper, this two-tiered strategy is implebe involved. The theoretical challenge is to quantitatively describe mented in a simpleway and applied the delocalizationof the transferring electron into and across the electronic interactions in a set of organic nonconjugated diene energeticallyunfavorable protein medium or the delocalization and dibenzosystems, a set of inorganicdiruthenium systems, and of a hole by an alternative mechanism? which intimately involves to electron transfer in cytochrome c. multiple exchanges among the electrons of the intervening Long-range electron transfer in complex environments can be medium. The former process is a form of tunneling in a complex theoretically treated in the very weak coupling, or nonadiabatic, medium,8 and the latter process is known as superex~hange.~J~ limit by a generalized version of Fermi’s Golden Rule,lI-l3 While the elementary theoretical foundations of indirect electron transfer are clear, it is evident that the challenge of k,, = (2a/h)lA12(FCWD) quantum biophysics is to perform realistic and quantitatively Thus the rate constant k, depends upon the product of two factors, useful calculationson a system of thousandsof electrons. Further lA)*and FCWD, arising, respectively, from the participation of motivation for a new look at this macromolecular quantum the electronic and nuclear motions in the kinetic process. The modeling task arises from the expectation that such tools developed Franck-Condon weighted density of final vibronic states, FCWD, for the biomolecular realm may become invaluablein the analysis is the quantum analogue of the activation factor in the original of complex molecular structures of synthetic origin as well. We Marcus rate expression, which describes and predicts the present a new level of macromolecular electronic structure analysis dependence of the rate upon the reorganization free energy and here, based on a well-defined strategy for the construction of a reaction e~ergonicity.1~ This factor is treated quantum mechanHamiltonian model, in which the macromolecular superexchange ically, semiclassically, or in the classical limit as appropriate to is treated as a consequence of the network of through-space and the degreeof vibrational excitationl”18 and need not be considered through-bond electronic couplings between sites within subunits

I. Introduction

QO22-3654f 9312097-5581%04.00/0 0 1993 American Chemical Society

5582 The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 further here. The superexchange electronic charge resonance energy A between the donor and acceptor electronicstates, on the other hand, is purely quantum mechanical. This is the critical distance and structure dependent variable that describes the transient electronic delocalization essential to electron transfer. In eq 1 the usual first-order electronic coupling matrix element, IHDA)~, where HDAis thedirect vacuum electroniccouplingbetween the zeroth order redox site wavefunctions, has been replaced with the infinite order superexchange coupling term, (AI2,for which an explicit expression can be derived.I1-I3 The formula for A is given in eq 3. Early (continuum) theories of long-range electron transfer were based on the WKB model for single-particle t ~ n n e l i n g ,which ~ ~ , ~treats ~ the intervening medium as a simple potential barrier, and predicts an exponentialdistance dependence of the direct or “through-space” donor-acceptor orbital coupling, falling off with the donor-acceptor distance, d, according to a constant, PDA,

where @DA is controlled by the ionization potential of the donor, EAA, at the transition state of the donor-acceptor system (at which moment it is also equal to the ionization potential of the reduced acceptor), and m is the mass of the electron. The problem is that the rate of falloff with distance given by eq 2 with any realistic nonphenomenological value of the barrier height is far more steep than that observed e~perimenta1ly.l~ Furthermore, eq 2, like any of the older simple exponential r ~ l e s , ’ ~ is J *an isotropicfalloff insensitive to the highly structured inhomogeneitys of the protein medium. Thus, to treat the electron transfer protein problem we turn away from the electron tunneling mechanism and turn instead to superexchange models. The required quantum mechanical method must be selected, or, in the present case, designed, with the topological features of the folded macromolecule in mind, particularly the importance of noncovalent interactions. Existing semiempirical methods are not designed, calibrated, or optimized to provide accurate descriptions of the delocalization of wavefunction tails nor the realistic representation of integrals at van der Waals contact range. In particular, with the atomic structure taken as given, a functional superexchange model will require shifting attention to longer length scales. Hence we have undertaken the assembly of a new model Hamiltonian. It is further important to appreciate that higher order perturbation theory, which is a natural approach to indirect or superexchange coupling, where HDA= 0, but where significant couplings HDI,, HI,]*, ...,HI,Ato and among the states of the intervening medium provide an electronic bridge from one redox site to the other, is most useful when the important bridge states are readily identified in advance,sv6but may be of restricted utility in treating electron transfer within the three-dimensional network of proteins. The distances between redox sites are often such that hundreds of atoms of peptide may be involved, with many noncovalent contacts between chains and a circuitous or nonexistent link along the covalent backbone; in such dense topologies simple combinatorials assures one that astronomically large numbers of paths may be r e l e ~ a n t .Most ~ importantly the identification of the nature of “bridges” is itself a key question in the biophysics of electron-transfer proteins, thus we wish to avoid restrictive assumptions at the outset, which are required by perturbation theory and few-path models. In the present nonperturbative or infinite order approach, the net charge resonance energy, A, arises as a full consequence of the threedimensional network of intertwined pathways of covalent and noncovalent superexchange coupling in the protein. The theory of A detailed below employsa schematic but carefully constructed matrix approach, which enables the unprecedented incorporation of the full three-dimensional multipath network involving the entire protein of cytochrome c in the nonperturbative evaluation

Gruschus and Kuki of A. While various implementations of this general strategy are possible, the approach detailed here is not an LCAO method but is rather based on an energetically inhomogeneous aperiodiclattice of sites located at each (nonhydrogen) atom and is referred to as IAL, the inhomogeneousaperiodiclattice model. An important measure of the success of this new approach is the ability to compute the value of the charge resonance energy A and not merely trends. The defining characteristics of this inhomogeneous aperiodic lattice strategy are the following: (i) The spatial energetic landscape of the hole superexchange medium is captured by an explicit calibration of each class of subunit comprising the macromoleculeso that diagonalizationof the subunit Hamiltonian subblock reproduces the experimental ionization potential (IP) of the functional group, thus avoiding a low accuracy recomputation of the already well-characterized experimental value. In the case of proteins, this means incorporation of as much as is experimentally known about IPS of the peptide backbone and side chains to calibrate the value of their diagonal Hamiltonian elements. (ii) The off-diagonal Hamiltonian matrix elements for both intra- and intersubunit interactions are calibrated by highaccuracy sources in the noncovalent and covalent limits. The noncovalent interactions are set by an interaction energy function that models the molecular orbital (MO) energy splittings of the u, r,and lone pair orbitals of small molecule pairs obtained by Gaussian 88 split valence calculations.21The covalent interaction, by which we mean the bonded-nearest-neighbor interaction, is set by calibration to experimentalsuperexchangeresults measured for smaller, fully covalently bridged molecules of well-defined geometry. (iii) The basis set is constructed solely for the purposes of representing subunits as a set of sites of interaction. This basis set may be defined at various levels of sophistication, but the entire macromolecule is included in a way allowingrepresentation of the entirevalence band. From the resulting Hamiltonian matrix the Green’s function required to compute A is then obtained by nonperturbative numerical matrix inversion. In the case of cytochrome c, this means inclusion of all 800+ heavy atoms in the cytochromes c in the matrix inversion algorithm with no restrictions or prior path search. Thus a distinguishing feature of this approach is that a schematic basis set can be employed, defined operationally in terms of the off-diagonal (eq 4) and diagonal elements (see Appendix). The IAL strategy does not compute subangstrom integrals or small-moleculeelectronicstructure but rather focuses on building up the macromolecular electronic couplings from the subunits, emphasizing realistic treatment of the crucial several angstrom interactions. The initial implementation given here of the general IAL strategy simply assigns every electron pair in the medium to the nearest heavy atom site. This gives rise to a strictly minimal basis set; that is, the number of basis functions used to express the Hamiltonian equals the number of occupied valence orbitals of the medium (half the number of valence electrons). Furthermore, anisotropicinteractions are ignored; a simple electronic interaction function is extracted according to (ii) above, but with the further simplification that the effects of orbital symmetry have been preaveraged. There are no explicit hydrogen sites, though their electrons are included in the total count of electron pairs. We expect that this relatively coarse-grained treatment of the details should describe the superexchange best in the longrange electron-transfer case, where the coupling arises from the sum of extremely large (combinatorial) numbers of paths. Neverthelesswe present the results from tests on relatively smaller molecules as well and emphasize that more refined implementations of the IAL strategy defined in (i)-(iii) above are clearly possible and are currently under development.

Model for Long-Range Electronic Superexchange Previous theories of long-range electron transfer have been developed in four major forms: path integral simulationmethods, high-order perturbation methods, pathway searching,and bridge wavefunction diagonalization methods. The first method requires a pseudopotential for the protein and has been applied to excess electron tunneling.8,22 This path integral approach requires multielectronfermionicexten~ions,*3~2~ which have the potential capacity to treat multielectron superexchange involving simultaneous hole (oxidized bridge) and excess electron states. The perturbation and diagonalization methods proceed through the explicit construction of a model Hamiltonian,whereas the pathway search, as Beratan and co-workers have developed it,25 evaluates potential pathways according to the number of links where each bond, hydrogen bond, or through-space link contributes a simple empirical attenuation factor. This homogeneous attenuation method, which treats all bonds between any two types of atoms in any compound as equivalent, provides useful qualitative mappings of tunneling pathways and has proven effective in giving ratios of electronic rate factors, A, agreeing with experimental res~lts.~ An~ interesting .~~ nonperturbativeformulationinvolving an iteration over all spatial sites has also been presented.28 Higher order perturbation methods, explored theoretically by Ulstrup and Kuznetsov,29JO were implemented on proteins by Ulstrup and co-workers3’ on selected sequences of amino acid fragments deemed a priori most important for electron transfer. Each sequence is a restricted linear perturbation chain in the sense that only a single term in the general perturbation sum is retained. Amino acid pathway calculations have also been performed by Broo & Larsson on azurin by bridge diagonalization with extended Hii~ke1.3~Marcus and co-workers included the effects of all pathways in a subvolumeof the interveningmedium by calculating, with extended Hiickel, the interactions between the donor and acceptor orbitals and fully diagonalized bridge molecular orbitals. Their noteworthy treatment of Ruthenium-modified Zn myoglobin is thus the first nonperturbative Hamiltonian-based electron-transfercalculation for a large domain within a protein, uses the same basis set and parameters as derived for smaller molecules and also provides reasonablevalues of the ratios of A.33 The present work is close in spirit to that of Siddarth and Marcus particularly in the emphasis on the simultaneous incorporation of very large numbers of pathways by full-matrix computations and in retaining the entire valence band, not just the frontier orbitals. The challengeremains, then, to calculate the charge resonance energy, A, and not simply a number whose magnitude only has meaning in relation to other numbers calculated by identical means. The required Hamiltonian must realistically incorporate the influence of the intervening medium upon long-range electron transfer, which implies the seemingly contradictory demands of a macromolecular computation at smaller molecule accuracy. The IAL Hamiltonian approach presented here sidesteps this difficulty of the monolithic approach by the use of the two-tiered strategyin which the small molecule and macromolecule electronic structure problems are separated. Note that all Hartree-Fock LCAO methods obtain their validation from the variational principle, yet it has been clearly established in many contexts that a rather extended basis set is required in order to obtain realistic results, particularly for interactions depending on wavefunction tails,2’934 and that minimal basis sets such as those one might contemplate using for a large system have little quantitative validity. Semiempirical methods such as extended Hiickel have a further difficulty for this application in that the ionization potentials can be several electronvolts in error relative to experiment. Both inaccuracies are avoided in the new twotiered strategy by loading both the diagonal and off-diagonal elementsof an effective Hamiltonianfrom higher accuracysources drawn from small molecule data.

The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 5583 The rest of the paper contains the details of the algorithm, followed by tests. Related superexchangecoupling elements, A, play critical roles in other long-range electronic phenomena for which experimental data are available. These include: (i) photoelectron spectroscopy measurementsof energy splittings of nearly degenerate double bond and aromatic orbitals coupled by nonconjugated bridge~,~5-’~ (ii) UV-vis absorption intensities of intervalence transitions in binuclear metal complexes held apart by rigid spacers,37(iii) remote heavy-atom fluorescence quenching,38and (iv) triplet-triplet energy t r a n ~ f e r . ~The ~ , latter ~ two processes are exchange-mediated, nonadiabatic, radiationless relaxation processes, analogous to electron transfer, while the first two are spectroscopic measurements,providing an important alternative source of superexchange data. Some of these molecular systems are just small enough to allow theoretical analysis using a variety of standard electronic wavefunction methods, but nevertheless we have proceeded to compute the IAL A in well-characterized experimental systems of types (i) and (ii). Specifically, we compute A for the polynorbornyl diene and dibenzo compounds of Paddon-Row and Jordad5J6and for the (Ru(NH3)5)2 dithiaspiro compounds of Taube and coworkers?’ where photoelectron spectroscopy and UV-vis absorption measurements have provided orbital splittinginformation on the bridge moleculeas well as intervalence transition intensities between the metal sites. These results test the performance of the IAL approach in the simpler covalent few-path frameworks and also serve to calibrate the covalent limit of the interaction function. With the required covalent parameter established, results for the electronic charge resonance energies for superexchangemediated electron transfer in three ruthenated cytochrome c proteins are then presented.

II. TheoryandMethods A. Nonadiabatic Rate Expressionfor ET in theSuperexchange Case. Define the initial and final states of the electron transfer (ET) reaction to be the product of the respective zeroth-order diabatic electronic state, $i,f, and the corresponding nuclear wavefunction. The diabaticelectronicstatesaredefined not with respect to the total Hamiltonian but rather with respect to the vibronically trapped electronic Hamiltonians, for initial and final electronic states in the absence of off-diagonal coupling between sites. In the initialdiabatic state the electronis localized on the donor, in the final state it is localized on the acceptor, and the net transition is made possible by the nonresonant transfer interactions between donor and protein and protein and acceptor via additional Hamiltonian terms not included in As described by Ulstrup, the total electronic Hamiltonian is equal to the diabatic Hamiltonians plus the coupling terms in the following way,3oH = Hi” 6 = H: Vr;though in the present long-range case, where the direct coupling ($dV&i) = 0, the distinction between Vi and Vf is inconsequential, and we may use a common V,which is dominated by off-diagonal couplings to and among the bridge. During the course of the nonadiabatic rate process, energy flows into and out of the reaction coordinate, and the diagonal energies of $i,f are modulated until the two are brought into resonance, whereupon a small probability of successful electron transfer occurs, as controlled by the square of the electronic charge resonance energy. Thus nuclear motions bring the initial system to the diabatically degenerate transition state, then, during the crossing event, the electronic interactions with and among the components of the protein matrix control the electronic delocalization, which gives rise to the successful electron transfer. A rigorous starting point for the quantum theory of long-range electron transfer in proteins is eq 1, originally developed to treat indirect processes in quantum scattering This equation and eq 3 provide the general nonperturbative superexchangerate expression for the transition probability per unit time in the

+

+

5584 The Journal of Physical Chemistry, Vol. 97, No. 21, 1993

nonadiabatic limit, where the result is valid to arbitrarily high order inHij, the protein-protein electronicinteractions(incontrast to approaching the adiabatic limit with high order in A4’). The electronic transition probability is proportional to JAJ2, where the key electronic charge resonance energy is found to

A = (+d~+,) + ( + ~ ( E -~ Ho + ie)-lq+i)

(3a)

(+dv+i)

(3b)

A = (qdqqi) i-x ( + d v + d ) ( E i o d

- Ed)-1

where Eio is the initial diabatic energy evaluated at the transition state and the positive infinitesmal E may be set to zero in the discrete matrix calculations here with no effect. In eq 3b, diagonalization of the Hamiltonian yields bridge eigenfunctions (+d], which may be inserted into the second term of eq 3a to give a conceptually useful formula for A explicitly in terms of a sum over couplings from +i to +d and from +d to +f. Either method of solution indicated by eq 3a and by eq 3b is very powerful as no approximate perturbation ~ u m m a t i o n and ~ ~ ,concomitant ~~ truncationgis required, though of the two numerical linear algebra problems, we choose the direct evaluation of eq 3a as more practical, since inversion of large matrices is considerably faster This inversion accomplishes a direct than diag~nalization.~~ numerical and nonperturbative evaluation of the elements of the Green’s functi0n,4~which describes the development of (A) I+i) (Eio - ZT)-14+l). The eigenvectors are thus not needed for the evaluation of A. In operation, the computation of the second term of eq 3a is accomplished by a vector-matrix-vector multiplication,where the vectors are the elementsof Vthat couple donor and acceptor to the bridge sites, and the matrix is (Eiol - H)-I where H is the matrix describing the bridge sites46and all interactions among them. Thus we further observe that the coupling from the same initial state to different target final states (+d can be quickly recomputed simply by changing the final vector (+$Vwithout recomputationof the Green’s function matrix (Eiol - H)-I, a potent advantage that we are exploiting in a variety of biomolecular applications. As protein reactions of interest are not self-exchangereactions; they have asymmetric environments that render it very difficult numerically to extract A as a small splitting between symmetric and antisymmetric eigenstates of the total Hamiltonian. There are no difficulties implementing eq 3 in asymmetric donor/ acceptor systems. While the splitting method is simpler, it is generally reserved for smaller, symmetrical systems.4s The Lbwdin partitioning method,49J0employed by Marcus and cow o r k e r ~and ~ ~by~Larsson,52 ~~ gives at first glance the identical result as eq 3, except that the Lijwdin partitioning method is strictly a time-independent electronic eigenvalue analysis, not a nonadiabatic rate theory, and hence the 6;’ appears there only as an approximation (first iteration step), rather than as a rigorous principle of electronic energy conservation as in the scattering expression here. High-order quantum scattering theory serves to provide a general time-dependent derivationof the nonadiabatic rate expression in the indirect coupling ~ a s e . ~ ~ -Aside l 3 from these formal distinctions in their derivation, they are the same nonperturbative matrix recipe. B. State Space for Superexchange and the IAL Hamiltonian. The superexchangeexpression eq 3 will be applied to three classes of processes: (1) photoelectron spectra (PES)splittings, (2) intervalence transfer, and (3) long-range nonadiabatic electron transfer. In the sections below we describe both the general algorithm and the specific implementations for the state space, the assemblyof the off-diagonal and diagonal Hamiltonian matrix elements, and the description of the redox or special sites for these three processes. In order to implement eq 3 to find A for the three transfer processes,one needs to describea model state space of intervening

+

-

Gruschus and Kuki medium, or “bridge”, states, with which the donor and acceptor states interact. Defining M to be the number of like-spin valence electrons in the bridge (M= 2500 in cytochrome c), each state involved is a configuration of M + 1like-spin electrons, including the one excess electron from the electron donor. Due to their indistinguishability, any of these electrons could be the one that ends up on the acceptor in the final state. Thus, as in any other configuration interaction problem, the states are properly Slater determinants, though as we will restrict the treatment to single excitations,thiscanbesimplified. Thebridgestateconfigurations are usually taken to be of two types, “hole” configurations characterized by occupied donor and acceptor sites, with one electron missing from a bridge valence orbital, and “electron” configurations, with no excess electrons at either the donor or acceptor sites but an excess electron in a bridge virtual orbital. Interactions between configurations of the same type have the possibilityof being appreciable,whereas matrix elementsbetween opposite types are usually negligible as this requires a higher order, two-electron process. Interactions of the initial and final configurations with both hole’ and electrons configurations can also be a p p r e ~ i a b l e .As ~ ~is well-known, the value of Eio (ED*) will control the relative dominance of the hole vs electron m e c h a n i s m ~ . ~In ~ Jour ~ calculations, we choose to include only the hole states, since, for the electron-transfer systems we have studied, the transferring electron energy is much closer to the top of the valence band than the bottom of the conduction band. For example, in the (Ru(NH3)& dithiaspiro compounds studied in this paper, the energy of the transfemng electron is approximately -6.3 eV (discussed below in section E), the highest bridge eigenvalue is -8.7 to -8.8 eV,54and the bottom of the conduction band for the bridge lies somewhere around 0 eV. In keeping with our relianceon experimentalnumbers for the bridge orbital energy levels, the 0 (* 1) eV value is taken from extensive experimental studies of quasi-free electron energies in fluid media.55 Multiple excitations are ignored, thus we adopt the view that the relevant superexchange configurations consist of the ground state of the system plus exactly one hole. The size of our Hamiltonian is then reduced from an astronomical number to M + 2, where again M is the number of like-spin valence electrons, or, equivalently, the number of bridge valence orbitals, and in this single excitation limit we can proceed to replace the more general configuration interaction approach by an effective interaction function within the one-particle orbital appro~imation.~~,5~ As with the extended Hiickel method, IAL depends upon appropriate calibrations without self-consistent field; but with the two-tiered strategy we go further in constructing a schematic lattice of interacting sites. The specific purpose of the IAL Hamiltonian is the realistic treatment of those site energetics and sitesite electronic interactions insofar as they affect long-range superexchange, where these matrix elements are computed (in the first tier) by prior analyses usingother basis sets, or experiment. The basis set employed in our current lattice model consists of assigning valence bond orbitals to the nearest heavy atom site, generating several valence bond orbitals (VBOs) per site, for a total of M VBOs. The mapping procedure first accounts for the heavy-atom/hydrogen VBOs by assigning one orbital to the heavy atom per attached hydrogen.. Next, one orbital per lone pair on the heavy atom is added. Finally, for every two bonds to other heavy atoms (double bonds count as two bonds), an orbital is assigned to the heavy atom. The total number of VBOs on the site is thus INT[no. of lone pairs + no. of attached hydrogens + (no. of bonds to non-H neighbors)/2]. The interaction between the VBOs is determined solely by the internucleardistancebetween the heavy atoms to which the orbitals are assigned, along with the parameters appropriate to the subunits to which they belong. The basis set is strictly minimal in that it spans only the valence band, and its diagonal energies are in the same range as the diagonal energies of valence bonds localized between two adjacent

Model for Long-Range Electronic Superexchange I

I

I

I

I

The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 5585 I

Y

Ab initio studies of diatomic systems have demonstrated that the use of highly flexible basis sets can provide a high-quality description of delocalization interactions across non-covalent distances as great as 13 A.34 Accordingly, the noncovalent behavior of the IAL off-diagonalelementsis set by an interaction energy function obtained by a fit to H F eigenvalue splittings of full multielectron Gaussian 88 ab initio calculations on small molecule pairs. Note that these computed orbital splittings do not arise from a simple molecular integral but rather incorporate the self-consistent interactions among the many electrons that constitute the molecular pair, including intermolecular delocalization and exchange. We have employed the triple-{ 6-31 1G (and some studies at 6-311++G) to compute the MO energy splittings arising from the delocalizationof VBO-like orbitals in pairs of ethane and pairs of formaldehyde molecules in various symmetrical configurations at up to 6-A intermolecular sepaI I I I I I u, rations. The comparison of the calibrated interaction function 0 1 2 3 4 5 6 7 with the computed splittings are given in Figure 1. The details nearest heavy atom distance (A) of the individual ab initio calculations, along with basis set Figure 1. An example IAL interaction function (dotted curve), depiction comparisons, have been given elsewhere.*’ of the covalent/noncovalent transition point determination, plotted along The calibration in the short-range, or covalent, limit is based with 6-31 1G basis Gaussian 88 calculations*’ of direct electronic on the calculation of the IAL superexchangecoupling elements delocalization interactions (one-half splittings) between ethanepair a-bond A for the smaller organic and binuclear molecules with wellorbitals (0)and between formaldehyde axial lone pair (a),*-bond (m), defined geometries. Various reasonable values for the throughand perpendicular lone pair (A)orbitals. Note that the ab initiocouplings bond coupling constant yo (see Figure 1 legend) are used; the were assembled onto a single distance abscissa by sorting as a function of distance between thenearest heavy atomsof each rcspective interacting comparison with the experimental A identifies the optimal molecular orbital. The choice of abscissa was the most effective in parameter yo, which is then used in the macromolecular superimposing the splitting curves from various 3D geometries2’ and applications of the IAL Hamiltonian. Thevariation of yocontrols defines theindependent variable dfor thecalibrationof the IAL interaction in part the valence band width and the ratio of covalent versus function V(d). The diagonal dashed line represents the noncovalent limit noncovalent interactions. Tables V and VI show the comparison of the interaction function, V,,(d) a - exp[-ld], where f is the average of A calculated by using values of yo from 1.0 to 3.0 eV with of the falloff constants of the interacting orbitals. The asymptotic falloff factorofagivenboundorbitaldecayingintovacuumisG = ( Z n ~ E j ~ / h ~ ) ~experiment /~ .60 if exchange effects are where EJ0 is the particular orbital The covalent limit yo is distinctly smaller than would beobtained diagonal energy and m is the mass of the electron. The intersection of by extrapolation back to one bond length of the asymptotic this dashed line, usin an average { = 2 A-l (values of 5 typically fall exponential from the ab initio noncovalent calibration, and hence with the through-bond coupling constant yo between 2.2 and 1.8 !-I), a transition is made from the exponential form valid at larger determines the transition point (tp). Thedotted lineis the IAL interaction distancesto the covalent yo coupling constant for distancesshorter function V(d),given in eq 4, for the case of yo = 1.0 eV. This interaction function is an exponential of a hyperbolafld), which serves to interpolate than a turnover distance a, between thecovalent and noncovalent asymptotes. acontrols the transition points’ and d (i0.5 A) modulates the sharpness of the covalent/ noncovalent transition.

atoms, rather than those of noninteractingatomicorbitals (AOs). The bonded nearest atomic neighbor AO/AO interactions, corresponding to the Hackel /3 parameters used in other treatmentsS* are thus subsumed into our VBO-type basis set diagonal energies. The through-bond interactions between neighboring VBO-type orbitals (for example, the matrix element between the two CC VBOs in propane) are bond-bond interaction matrix elements, referred to here as yo, and cause a single broad band to form, which is the valence band. While it is not required in the general IAL strategy to employ such a schematic basis to represent the sites of interaction in each subunit, the present construction makes the method readily applicable to quantum modeling of large macromolecular systems. C. The InteractionFunctionfor Off-Diagd Matrix Elements. The IAL strategy focuses on the two properties deemed most essential for long-range superexchange prOCeSses, the site-site, off-diagonalinteractions that lead to delocalizationof the redox site wavefunction tail (this section) and the variations in local energetics of the superexchangemedium that guide and modulate this delocalization, determined for the most part by the orbital diagonal energies (section D). The magnitude of the off-diagonal electronic delocalizationinteraction energies is modeled as a single function obtained by calibrating the interactions separately in the covalent and noncovalent limits. The covalent limit in the IAL model is the interaction between nearest bonded neighbor sites and is controlled by an asymptotic short-range coupling parameter,yo,while the asymptotic noncovalent exponential slope is set by standard electronic structure methods as follows.

This largely exponential V(d) expression is the IAL interaction function for off-diagonalmatrix elements. Longer-range,direct interactions are complicated by the possibility that a third, occupied orbital has come between the two interacting orbitals, and so we have chosen a cutoff for direct interactions at six angstroms. Theeffectofthecutoffwasexamined;seethefootnote to Table VII. Further details of the parameters {, 2, and 6 of the interaction function for the total distance range d C 6 A, shown in Figure 1, are given in the figure legend. At short range, the aperiodic lattice of sites can only approximately represent real interacting valence bonds, which are directed and bond-centered. Bond-centered VBOs would be shifted by half a bond from the present simpler atomic site-centered VBOs.6l VBOs assigned to the same heavy atom often originated from VBOs with different symmetries; these VBOs are simply assumed to have no interaction. This “same-atom-no-interaction”rule might be said to reflect the obviouseffect of symmetry on interactions in conjugated systems between u and T bonds, or between lone pairs and T bonds, though it is less appropriate for two Q VBOs mapped to

Gruschus and Kuki

5586 The Journal of Physical Chemistry, Vol. 97, No. 21, 1993

TABLE I: Pol orbornyl Bridge Parametere with Top aod Bottom VIlleace Band Energies for Different Values of the Through-Bond &pling Constant 70in the IAL Interaction Energy Function C diagonal energy,‘ eV -17.69, yo = 3.0 eV -14.08, yo = 1.5 eV -15.17, yo = 2.0 eV -12.97, yo = 1.0 eV eigenvalues,eV 1st last 1st last 1st last 1st last diene (bridge only) n=O n=l n=2 n=3

-9.97 -9.85 (-9.65) -9.59

-22.99 -26.04 -27.58 -28.39

-9.98 -9.70 (-9.66) -9.68

-28.12 -32.04 -33.97 -34.97

-9.84 -9.70 (-9.6 5) -9.53

-32.32 -36.68 -38.76 -39.82

-9.99 -9.82 (-9.65) -9.55

-39.59 -44.37 -46.55 -47.65

-9.96 -9.86

-22.97 -26.02

-9.98 -9.70

-28.10 -32.01

-9.84 -9.70

-32.28 -36.64

-9.98 -9.82

-39.54 -44.32

dibenzo (bridge only)

n=O n=l

The orbitals assigned to the carbons of the n = 2 diene compound (see Figure 2), minus the four “ene” carbons, were used to generate an IAL Hamiltonian. The diagonal energy of this Hamiltonian was then uniformly adjusted so that, upon diagonalization, the highest eigenvalue equaled -9.65 eV (in parentheses above), which is the first peak in the PE spectrum beneath the diene peaks.35 Then = 2 compound was chosen because examination of its PE spectra showed the most readily identifiablebridge IP of all the compounds. While the IAL strategy normally calls for calibrating the individual subunits of the medium, in this case the norbornane spacers of the bridge, the fact that the subunits in the n = 2 compound are identical means the same parameters will be used for them anyway. Ideally, calibration of then = 1 compound to its bridge IPshould give the same diagonal orbital energy. As per the orbital assignment rules, the methylene carbons each received three orbitals, two for their attached hydrogens plus a total of one orbital from their two bonds to neighboring heavy atoms. The methine carbons each received two orbitals (rounded down from 2.5) for their one attached hydrogen plus a total of 1.5 orbitals from their three bonds to neighboring heavy atoms.

Figure 2. Structures of’ the polynorbornyl compounds. For the diene compounds n = 0, 1 , 2 , 3 and for the dibenzo compounds n = 0,l. The atomiccoordinates are from the MNDO optimized geometry;the MNDO calculationswere carried out by the authors on all the diene and dibenzo compounds. X-ray crystal structures are also available for the diene compounds.35

the same atom. In addition, the interaction function approach constitutes an s-orbital-likemodel, though an important distinction is that these simplifyingrules are not dependedupon in computing the subunit energetics, which are fixed instead by calibration via ionization potentials (section D). As is clear from Figure 1, the interaction curve describes an average over various orbital symmetries. D. Calibration of the Bridge Orbital Diagonal Energies. The specific diagonal energies of each type of functional group constituting the medium are calibrated to functional group ionization energies, not to atomic ionization energies. The diagonal site energies within a group were uniformly adjusted so that upon diagonalization of the subblock of the IAL Hamiltonian corresponding to the group, the highest energy eigenvaluewould match the experimental IP of that functional group.62 The algorithm for the refinement of the diagonal VBO energies is detailed in the Appendix, with the associated Figure 5 . (1) Calibration of the Orbital Diagonal Energies of the Polynorbornyl Compounds. For the calculation of the superexchange coupling elements A, the superexchange medium, or bridge, consists of the atoms of the particular polynorbornyl molecule minus the alkene and arene subunits whose interaction is being studied. The results and specific information on the calibration of the diagonal energies are given in Table I (see also Figure 2). (2) Calibration of the Orbital Diagonal Energies of the Wthiaspiro Compounds. The study of the dithiaspiro compounds

involves two sets of compounds; for the sulfur lone pair PES splittings in the disulfides, the bridge consists of the intervening cyclobutyl rings, and for the intervalence transitions in the (Ru(NH3)& binuclear complexes, the bridge is the same except that the sulfurs are included. Table I1 (see also Figure 3) contains the results and specific information regarding the calibrations for both cases. (3) Calibration of the Orbital Diagonal Energies of Protein Functional Groups. The cytochrome c protein superexchange medium provides by far the greatest variety of functional groups of the systems studied. These functional groups have been individuallyconsideredin three categories: easilyionizablegroups, non-easily-ionizablegroups,and peptide backboneor “bulk”.The first category consists of the aromatic-, carboxylate-, and sulfurcontaining side chains (Trp, Tyr, His, Phe, Cys, Met, Asp, Glu), while thesecond is made upof the nonpolaraliphaticand positively charged side chains (Ala, Leu, Val, Ile, Pro, Arg, L Y S ) . The ~~ remainder of the protein, the amide backbone and polar but nonionicaliphatic sidechains (Asn, Gln, Ser, Thr, Gly),constitutes the third category. The orbital diagonal energy calibration of each easily ionizable groups follows a two-step process similar to the dithiaspiro bridge calibrations including the sulfurs (see Table 11). For each type of side chain, a selection of residues of that type is extracted from the coordinate data for the full protein. The alanine fragment is first calibrated with a separate (deeper) diagonal value from that of the specific functional group in the side chain. Then with the side chain re-attached, the side chain diagonal elements are calibrated. Figure 4 depicts the two-step process for tryptophan; see also the Appendix. Specifically, the easily ionizable group atoms are those of the indole, phenol, benzene, and imidazole rings for the aromatic residues, the sulfur atoms of methionine and cysteine, and the carboxylate oxygens of the carboxylate side chains, including those of the propionate groups of the heme. The remaining carbonsof thesulfur- and carboxylate-containingamino acids are assigned the alanine orbital diagonal energy in the calibration. The calibration of the non-easily-ionizable groups closely parallels that of the easily ionizable groups. First the orbital diagonal energies of the glycine portion of the residue are calibrated and held fixed. The remaining deeper orbital diagonal energies of all the atoms in the aliphatic or positively charged side chain are then calibrated to reproduce upon diagonalization the a-bond MO onset in the PES, which lies deeper than the amide IP. Table I11 summarizes the results for the protein functional group calibrations with yo = 1.5 eV. All the numbers given are

The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 5587

Model for Long-Range Electronic Superexchange

TABLE II: Comparison of Dithiaspiro Bridge Eigenvalues for Different Values of the Through-BondCoupling Constant 70in the IAL Interaction Energy Function with Photoelectron Spectra Peak Values for Dithiaspirane Disulfides ~

~~

~

yo = 1.0eV

eigenvalues (ev) 1st compd I -8.84 compd I1 (-8.75) compd I11 -8.76

2nd -9.48 -9.34 -9.20

~~

~~

yo = 1.5 eV

3rd -11.49 -11.35 -11.09

1st -8.84 (-8.75) -8.79

2nd 10.74 -10.20 -9.67

yo = 2.0 eV

3rd -11.79 -11.06 -10.90

1st -8.54 (-8.75) -8.96

2nd -10.93 -9.89 -9.37

~

yo = 3.0 eV

3rd -11.85 -11.47 -11.23

1st -8.06 (-8.75) -8.61

2nd -9.73 -8.78 -8.92

3rd -13.15 -11.25 -10.79

experimental PE spectra peaks“ 1st peak 8.71 8.75 8.75

2nd peak 10.04 9.55 9.45

3rd peak =ll =ll =11

bridge diagonal orbital energies, eV carbonb sulfurc

yo = 1.0eV

yo = 1.5 eV

yo = 2.0 eV

yo = 3.0 eV

-14.14 -1 1.49

-1 5.29 -1 3.86

-16.75 -14.80

-19.42 -1 6.45

Values taken from ref 54 for the disulfide compounds, where the third “peak” actually corresponds to the initial leveling-off of the shoulder at the onset of the band of lower energy levels. Following convention, positive peak values are reported, though their negative should be compared to the eigenvalue energies. The orbitals assigned to the carbons of compound I1 (see Figure 3), without the sulfurs, were used to generate the IAL Hamiltonian. Two orbitals were assigned to the quaternary carbons and three orbitals to the methylene carbons, as per the orbital assignment rules. The diagonal energy of this Hamiltonian was then uniformly adjusted so that, upon diagonalization, the highest eigenvalue equalled -1 1 eV, the approximate onset of the nonsulfur MO energies in the disulfide ~ o m p o u n d . ~Using ~ , ~ ~the diagonal value predetermined for the carbons, the IAL Hamiltonian was generated by using all the bridge orbitals of compound I1 (including the three orbitals assigned to each of the disulfide sulfurs). The diagonal energy of the sulfur orbitals was then adjusted so that the highest eigenvalue of the diagonalized Hamiltonian matched the experimental IP of disulfide bridge compound 11, 8.75 eV (this eigenvalue is in parentheses). The corresponding disulfoxide IP is effectively the same, 8.78 eV. CONDUCTION BAND = 0 eV

= 0 eV

VALENCE

Figure 3. Structures of the dithiaspiro bridges and Ru(NH3)s ligand conformations. The atomic coordinates for the bridges are from the MNDO optimized geometries, calculated by the authors. For the disulfide PES splittingstudy, the MNDO geometries werecomputed on thedisulfide compounds. In the case of the calculations of the binuclear intervalence transfer superexchange coupling, the MNDO geometry calculations were performed for the corresponding dithiaspirane disulfoxides. Two rotational conformationsof the Ru(NH3)5ligands were used in the intervalence transition calculations: (1) the eclipsed conformation where three of the nitrogens, the ruthenium, and the sulfur are in a plane bisecting the sulfur-containing ring and (2) the staggered conformation, rotated 45’ around the R u S bond from the eclipsed conformation. The ligands on both metals were given the same rotational conformation. The Ru-N distance of 2.10 8, was taken from the X-ray crystal s t r ~ c t u r e , while ~’ the R u S distance of 2.19 A was taken from ref 37b. The pyramidal ligation of the ruthenium was accomplished by extending the S-0 bond length to that of R u S .

the averages of all the particular residues studied; note also footnote f on the averaging for non-easily-ionizable groups. E. Donor and AcceptorOrbitals. The special terminal orbitals $i,f are modeled as the set of sites located on the heavy atoms of the groups and bearing a uniform distribution of spherical VBOs, one per site. Normalization constants are calculated for each of these $ i , ~ MOs, including the correction for orbital overlap, obtained with the simple overlap formula for Slater-type 1s orbitals. The interaction element between the donor or acceptor orbital and a given bridge orbital is calculated by finding the interaction of each donor or acceptor atom with the bridge orbital by using the interaction function, summing over atoms, and multiplying by the normalization factor. The introduction of accurate experimental energies for the redox sites requires an accurate source of data, which we take to be experimental reduction potentials. These are readily available and all manner of detailed variations due to environmental or ligation effects are straightforward to assess by electrochemistry. In order to translate these data into the IP values required, we use the following empirical correlation

-7.90 eV

Figure 4. The two-step orbital diagonal energy calibration process for easily ionizable protein residues using tryptophan for an example (also see Appendix). The first step truncates the residue side chain at CS, leaving the alanine fragment, and uniformly adjusts the value of the orbital diagonal energy elements, Hii, for the alanine fragment, until the highest eigenvalue of the diagonalized Hamiltonian matches the experimental alanine IP of 9.85 eV, which requires an orbital diagonal energy of -15.2 eV. With the Hii for the VBOs in the alanine fragment held fixed at this value, the IAL Hamiltonian for the entire residue is then constructed. In the second step only the Hii elements of the VBOs assigned to the easily ionizable group atoms, the indole atoms in the tryptophan case, are refined. The higher diagonal energy of -12.8 eV is found to yield the required match between the highest subunit eigenvalue and the tryptophan IP of 7.9 eV. Steps 1 and 2 are shown on the left and right, respectively, as the development of the corresponding bands from the optimized diagonal VBO Hii inputs in the center. Due to the strictly minimal basis set (see section B) the IAL Hamiltonian only produces eigenvalues corresponding to occupied MOs; no conduction band eigenvalues are produced (shown as empty above). The off-diagonal interactions (defined by eq 4), which develop the valence band, represent averaged direct electronic couplings between occupied valence bond orbitals, from triple-{ Hartree-Fock model calculations.

E,/,(NHE) = 0.92IP - 5.66

(5a)

IP = (E,/,(NHE) + 5.66)/0.92 where the correlation is between vertical IP and the thermodynamic (adiabatic) single-electron EI,2(NHE);70see Table IV. Solvation effectsand the differencebetween vertical and adiabatic IP values are averaged in such a correlation. As a separate check, the ferrocene E1,,(NHE) = +OS9 to 0.64 V in acetonitrile yields IP = 6.79 to 6.85 eV, which is not exact but quite reasonable (the IP is 6.91 eV; see also ref 77). (1) Smaller Compounds. For the polynorbornyl compounds and the dithiaspirane disulfidecompounds, the atoms of the donor

Gruschus and Kuki

5588 The Journal of Physical Chemistry, Vol. 97, No. 21, 1993

TABLE III: Orbital Diagonal Energies for Protein Functional Group (yo = 1.5 eV) easily-ionizable groups

residue

orbital diagonal energy, eV

TrP TYr Cysd Met Phe His average

-12.8 -13.7 -12.2 -12.5 -13.5 -13.9 -13.1'

non-easily-ionizable groups

IP, eV

residue

orbital diagonal energy, eV

7.9" 8.6" 8.6" 8.6" 9.1" 9.1'

Ala Val Pro Ile Leu

-17.8 -17.2 -18.6 -16.9 -18.6

peptide backbone

U-MO onset, eV -12.6" -12.lb -12.0" -11.9' -11.8'

-17.W

residue

orbital diagonal energy, eV

IP," eV

GlY Ala

-15.7 -15.2

10.0 9.85

-15.48

" Reference 65.

Reference 66. The IP for histidine, 9.0 eV,76was not known by the authors at the time of the ~ l i b r a t i o n .The ~ ~ value employed in the calibration is adapted from pyrrole.68 The only cysteines in the cytochrome c proteins studied are covalently linked to the heme at the vinyl group positions. These two additional carbons are included in the calibration. e As there is no experimental information on the IPSof negatively charged side chains, the carboxylate oxygens are assigned the average of the easily-ionizable group orbital diagonal energies. f As there is no experimental information on the IPS of positively charged side chains, the positively charged side chams are assigned the average of the non-easily-ionizable group orbital diagonal energies. The IP of H2O is 12.6 eV, so trapped water orbitals are also given the average energy. Furthermore, due to the relatively less accurate values of experimental u-bond MO energy band onsets and the less precise nature of the hydrophobic aliphatic calibration procedure, all hydrophobic aliphatic groups have been assigned the average energy. All remaining orbitals (backbone and hydrophilic aliphatic side chain orbitals) have been assigned the average peptide backbone value. Note that the actual average values calculated for alanine and glycine were -1 5.70 and -15.16 eV, so their aierage is rounded to -15.4 eV.

TABLE IV: Ionization Potentials from Redox Potentials of Electron-Transfer Sites complex RU(NH3)63+/2C

Ru(NH&-His3+I2+ Zn-porphyrin'+/O Zn-porphyrin'+/*

+

electrolyte

NHE,V

IP,"eV

0.3 M Na[CH3SO3] 0.1 M Na[BF4] 1 M KCI 0.036 M Na[C104] 0.1 M Na[C104]

+0.057 +0.051 +0.095 +0.214 +0.080 1.09 -0.62

6.21 6.21 6.26 6.38 6.24 7.34 5.48

+

IP = ( E , I ~ ( N H E ) 5.66)/0.92adaptedfromEl/2(Ag/Ag+) = 0.92IP -6.20 byusingEl/2(NHE) zE1/2(Ag/Ag+) +0.54V,ref70. Reference 69. Reference 2b. Reference 2c.

'

and acceptor orbitals are the carbons of the alkenes and benzenes and the sulfurs of the disulfides with equal amplitude for each atom. For the Ru(NH3)5 ligands, the ruthenium and the five nitrogens were given amplitudes taken from molecular orbital populations of a calculation done on R u I I ( N H ~ ) ~ . ~ ~ The energiesEio for the donor and acceptor orbitals were taken from the experimental IPS of appropriate compounds when available. For the alkene orbitals the IP of cis-Zbutene, 9.13 eV,68 was used, and for the benzene orbitals the IP of o-xylene, 8.555 eV,6*was used. For the sulfurs of the dithiaspiranes, the IP of dimethyl sulfide, 8.68 eV,68was used. The energy for R u I I ( N H ~ was ) ~ derived by interpreting several reduction potentials in different solvents through the empirical formula (eq 5 and Table IV). The error of the empirical formula appears to be h0.3 eV, hence two more values of Eo = 0.3 eV above and below the average derived IP of 6.3 eV were used in additional calculations from which the sensitivity of the superexchange coupling to Eio could be determined. (2) Cytochrome c. The two electron transfer sites in the cytochrome c derivatives under investigation are the interior zinc porphyrin and a pentaammine Ru complex coordinated to a surface histidine. The zinc porphyrin plays the role of donor in the forward transfer upon photoexcitation and nanosecond relaxation to the operative porphyrin triplet state. Its *-radical cation becomes the acceptor in the thermal back transfer. The atoms of the zinc porphyrin orbital include the zinc, the four pyrroles, and the four meso carbons. The methyl, alkylcysteine, and propionate side chains are considered instead part of the superexchange medium. The atoms of the pentaammine Ru complex include the five nitrogens and ruthenium of the complex and the five atoms of the coordinated imidazole ring. All atoms were assigned equal amplitude, and the normalization constants and interactions with bridge states are determined as described above.

The reduction potentials of the ruthenium and porphyrin sites are used to determine the relevant donor and acceptor orbital energies as described above for the smaller moleculesand detailed in Table IV. In a given electron-transfer system, the vertical redox potential of the donor/acceptor system at the essentially degenerate transition-state configuration defines the energy of the active initial and final electronic orbitals. We will refer to this transition-state electronic energy Eio as EDA,the energy of the transferring electron at the moment of degeneracy between donor and acceptor orbitals. The donor electron energy and the acceptor hole energy, determined through eq 5, and the reorganization energies around the respective sites together determine this redox potential at the transition-state configuration (see footnote 44 in ref 9). With the simplifying assumption of equal contributions from the two sites to the reorganization energy, the forward and backward transfer electronic energiesEDAare taken as the arithmetic average of the Ru complex and appropriate zinc porphyrin orbital energies, giving -5.9 eV for the excited-state reaction and -6.8 eV for the ground-state reaction. This systematic procedure for the estimation of EDAmay be applied to any condensed-phase electron transfer in which the reduction potentials are known. Tests of the sensitivity of the result to the precisevalue of EDAhave also been carried out as described below. 111. Results

A. Orbital Splittings. (1) Polynorbornyl Diene and Dibenzo Compounds. The superexchange coupling elements A between alkene *-bond orbitals and between benzene *-system molecular orbitals were calculated by inverting the IAL Hamiltonian and the use of eq 3a and compared with one-half the experimental orbital energy splitting from the PES. These molecules position the *-orbitals apart with zero to three intervening norbornane units for the diene molecules and zero and one norbornane unit for the dibenzo molecules (see Figure 2). The orientation of the alkene and benzo groups was designed to enhance noncovalent interactions with the norbornanespacers, and the resulting pattern of superexchange geometry is known as laticyclic hyperconjug a t i ~ n This . ~ ~ combination of important noncovalent interactions alongside the covalent interlinkages provides a valuable system for testing the IAL Hamiltonian and calibrating its throughbond coupling constant. Table V presents the IAL Hamiltonian results for the superexchange coupling eIements A, calculated by using four values of the through-bond coupling constant yo. Good overall agreement with experiment is achieved for yo equal to 1.5 and 2.0 eV, with most values of A slightly less than the experimental half-energy splittings. It is striking that, despite a 3-fold increase

Model for Long-Range Electronic Superexchange

The Journal of Physical Chemistry, Vol. 97, No. 21. I993 5589

Comparison of IAL Superexchange Cou Elements A to MNDO and Experimental Orbital Energy Splittings in DithiasDirane Disulfides and Polvnorbornvl Diene an DI lizo Comoounds yo = 1.0eV yo = 1.5 eV yo = 2.0 eV yo = 3.0 eV HDA" A HDA" A HDA" A HDA' A MNDOd I/&!? exptl I / @ dist, A

TABLE V

Pie

-0.1OC compound I -0.038 compound I1 -1.4 X 10-' +0.10 -0.080 compound 111 -4.9 X

Sulfur Couplings,beV -0.5lC -0.046 -0.31e -0.049 -0.4SC -0.054 -1.7 X l t 3 +0.28 -1.8 X l t 3 +0.27 -2.0 X lP3 +0.033 -5.8 X -0.20 -6.2 x 10-5 -0.11 -6.9 x 10-5 -0.083

n=l n=2 n=3

-0.59 -0.39r +0.16 -5.7 X -5.6 X 1Cr5 -0.055 +OB067 -8.4 X

-0.73 -0.48c -6.8 X lC3 +0.28 -0.14 -6.7 X -1.0 X 10-6 +0.034

n=O n=l

-0).27c -0.33 -3.4 x l t 3 +0.044

-0.40 -4.1 X

n=O

~~

Diene Couplings, eV -0.80 -0.48c -7.2 X lo-) +0.24 -7.1 x 10-5 -0.12 -1.1 X 10-6 +0.050

-0.90 -7.9 X 1 C 3 -7.8 x 10-5 -1.2 X 10-6

fl.2 x 1 w i 2 . 3 x 10-3 f 7 . 1 X lo4

d~0.67~ f0.40C f0.35e

-0.53c +O. 16 -0.067 +0..0028

-0.32 -0.066 -0.021 -0,0072

-0.63f -0.26f

1O.OY

-0.3SC +0.042

-0.14 -0,027

-0.5or -0.19

-0.19

4.88 7.10 9.32 2.99 6.16 9.18 11.90

Dibenzo Couplings, eV -0.32c -0.43 -0.34c -0.49 +0.066 -4.4 x 10-3 +0.063 -4.9 x

10-3

3.189 6.238

a Direct vacuum coupling between the donor and acceptor, using the IAL interaction function ignoring the 6-A cut-off. * The IAL A was calculated with the disulfide MNDO optimized geometry and eq 3a. The direct vacuum term is included in A (interacting orbital distance 1 6 A). MNDO ~.~~ values taken from photoelectron spectra calculations done by the authors. e Experimental values taken from photoelectron ~ p e c t r a . )IExperimental with the sign inferred from ab-initio calculation^.^^^^^ 8 Nearest benzene-benzene carbon atom distance.

TABLE VI: Superexchange Intervalence Transition Coupling A for Dithiaspiro (Ru(NHJ)& Compounds for Different Couplings Constants yo in the IAL Orbital Interaction Function direct vacuum HDA,"cm-I coupling IAL A$ cm-' exptlC ext Huckel distances: A, stag. eclipsed A, cm-I A, cm-I stag., eclipsed stag. eclipsed const YO, eV Compound I 1.o 1.5 2.0

-54.5 -67.6 -74.4

1.o 1.5 2.0

-1.53 -1.89 -2.08

-1.63 -2.02 -2.22

1.o 1.5 2.0

-0.258 -0.320 -0.352

-0.284 -0.352 -0.387

-50.9 -63.1 -63.4

+141 (26) +114 (1 1) -222 (141)

+132 (20) +76 (23) -325 (167)

f207

171

7.62 4.78, 5.40

178

57

10.09 8.53,8.16

137

20

11.92 9.15, 8.98

Compound I1 -28 (4) +11 (23) +191 (80)

-29 (3) +24 (30) +227 (99)

Compound I11 +2.8 (0.6) -22 (16) -90 (39)

+2.2 (1.4) -29 (19) -100 (42)

The direct vacuum coupling HDAis not included in the IAL results. The IAL calculations were done with the Ru(NH3)s energy equal to 6.3 eV (see Table IV). Calculationswere also done at 6.0 and 6.6 eV, and the maximum deviations from the 6.3-eV results are given in parentheses as a measure of the sensitivity of the IAL results to the Ru(NH3)s energy. Reference 37b, but with the A values corrected by using the Ru-Ru distances reported in the last column above which arose from the MNDO optimized geometries for the bridges (see Figure 3). The first value is the Ru-Ru distance. The two distances beneath are the nearest nitrogen-nitrogen distances in the staggered and eclipsed conformations. in yo, the IAL results do not respond monotonically but rather, for n # 0, peak at yo equal to 1.5 and 2.0 eV. Note that the IAL bridge diagonal matrix elements (section IID) respond to changes in eq 4 and that, while the highest bridge valence eigenvalue is held constant by the recalibration procedure, the average bridge eigenvalue is not. The nonmonotonic behavior of A can be understood to result from the valence band broadening obtained from the corresponding yo = 3.0 eV bridge calibration and the energy denominator in eq 3b. The broadening results in lower average bridge energies as reflected in both the lower IAL diagonal matrix element of the bridge orbitals and in the deepest valence eigenvalue of the subsequent diagonalized bridge Hamiltonian ("C diagonal energy" and last eigenvalue in Table I, respectively). The signs of the IAL superexchange elements A indicate the presence of a nodal structure in the superexchange within the medium, even though the lattice model used here employs basis functionswith spherical symmetry.@ Disagreementof these signs with those from more sophisticated MNDO and ab initio STO3G calculations indicates that an accurate determination of the sign requires inclusion of appropriately directed p-type orbitals. This is a limitation of the present coarse-grained no-hydrogens implementation of the IAL strategy or any quantum model that stops short of 'atomic resolution". For quantitative comparison with a semiempirical SCF-MO alternative, one-half the orbital energy splittings from the MNDO Hamiltonian were also computed and are included in Table V.

(2) DithiaspiraneDisulfides. The IAL superexchangecoupling elements calculated between sulfur lone pair orbitals in dithiaspirane disulfide compounds with two to four rings (Figure 3) are also given in Table V and compared with one-half the PES ~ p l i t t i n g s .In ~ ~Table V the IAL results for A are presented for the four values of yo, and, in addition, Table I1 contains IAL eigenvalueswhich can be directly compared with the experimental PES peaks. (In symmetrical compounds with larger couplings ( eV), the direct eigenvalue results (Table 11) suffice, and A calculation by eq 3a (Table V) is not required.) While the underestimation of A for the sulfur superexchange coupling elements is on the order of a factor of 2, the splittings from Table I1 tend to slightly overestimate the experimental splittings, and, overall, there is still reasonable agreement for yo = 1.5 and yo = 2.0 eV. Also, the same nonmonotonicbehavior of A as a function of the parameter yois seen. The results from MNDO calculations are presented for comparisonand can be seen to be in considerable error, severely underestimating the superexchange coupling. B. Superexchange Coupling for Intervalence Transitions in Ru11(NH& and Ru1I1(NH&Dithiaspiro Compounds. The superexchange coupling elements A between Ru(NH3)s sites in the coordinated dithiaspiro compounds (see Figure 3) were calculated by the IAL Hamiltonian and compared with values of A extracted from intervalencetransition peak intensities.3'b Calculationsdone by using three values of yo are presented in Table VI. Fixing both ligands in either the eclipsed and staggered conformations

5590

Gruschus and Kuki

The Journal of Physical Chemistry, Vol. 97, No. 21, 1993

TABLE VIk Charge Resonance Energies A for Electron Tunneling in Rutbennted Cytochromes c (YO = 1.5 eV) redox site edge-edge dist, A His39" Hi~33~ Hi~62~

13.0 13.1 15.6

EDA= -5.9 eve direct vacuum IAL HDA,cm-I A, cm-I

EDA= -6.8 eVf direct vacuum IAL HDA,cm-I A, cm-l

exptld A, cm-I

-2.3 x 10-3 -4.1x 10-3 -1.4 X 1 W

-8.9 X 10-4 -1.5 X 1 C 3 -4.1 x 10-5

h0.21 *0.12 *0.012

+0.20 +O. 18 -0.018

-0.74 +0.28 +0.019

Yeast Candida krusei cytochrome c derivative.72b Horse heart cytochrome c derivative.72b Triple mutant cytochrome c derivative expressed via yeast Saccharomyces c e r e ~ i s i a e . ~ Reference *~ 27. Calculations also done for EDA= -5.8 eV give sensitivities of A to EDAof +0.13 cm-I/eV, -0.17 cm-l/eV, and +0.007 cm-'/eV for the His39, His33, and His62 derivatives. (Calculations also done for EDA= -6.9 eV give sensitivities of A to EDA of +5.2 cm-l/eV, +0.27 cm-'/eV, and -0.25 cm-I/eV for the His39, His33, and His62 derivatives. Calculations were also redone by using a 7-A cutoff, which gave the same diagonal parameters as the calibrations done with the 6-A cut-off. IAL results for the 7-A cut-off and E D A= -5.9 eV are 0.20 cm-I, 0.18 cm-I, and -0.019 cm-I, and for EDA= -6.8 eV, -0.72 cm-I, 0.32 cm-I, and 0.030 cm-I for the His39, His33, and His62 cytochrome c derivatives, respectively. Note practically no effect occurred for EDA= -5.9 eV, while for the longest distance case a modest change can be seen for E D A= -6.8 eV.

enabled an examination of the impact of varying donor-acceptor orientations in a simple system with a rigid bridge compound. Additional calculations using ligand energies Eio f 0 . 3 eV from the average deduced RuI1(NH3)5 IP (e6.3 eV, see Table IV) provided a measure of the sensitivity of the results to EiO. The sensitivity of the results to Eio is quite substantial; for instance, for compound I1 with yo = 1.5 eV, the change in A with a 0.3-eV change in Eio is greater than the magnitude of A, and the sign of A does indeed change. This sensitivity, especially in this particular case, is the inevitable consequence of the nodal structure of the delocalized hole state. The variation of the superexchange coupling with respect to the donor-acceptor rotational conformation is less drastic than for changes in EiO, though still significant. No choice of yo is uniquely superior in producing the best values of A. The choices yo = 1.5 eV and yo = 2.0 eV provide calculated values of A that bracket the experimental data. On the other hand, the results for yo = 1.0 eV fall more than a factor of 10short for compound 111, which indicates insufficientcoupling. The overall performance in Table VI is not as good as in the case of the organic compounds, probably reflecting the greater need to focus on the detailed redox site wavefunctions. Comparing the accuracies of the IAL A for all the intervalence transition and orbital splitting results above, the value of yo = 1.5 eV was selected for use in the subsequent macromolecular calculations. Given the overall consistent performance of yo = 1.5, and accepting a factor of 2-3 accuracy from the inhomogeneous aperiodic lattice model, a higher precision calibration of yo is not justified at this level of IAL implementation. Extended Hiickel results are shown for comparison. Their agreement with experiment is excellent. It is curious that this semiempirical method does so well here, while the more sophisticated MNDO method did so poorly in predicting orbital splittings (Table V). C. SuperexchangeCoupling Elements for Electron Transfer in Cytochrome c. Three different protein structures were used in our study, one horse heart cytochrome c and two yeast cytochromes c. The electron-transfer rates for surface histidine-ligated Ru derivatives of each of these proteins have been measured e ~ p e r i m e n t a l l y .The ~ ~ ~three , ~ structures have ruthenated surface histidines at residue 33 for the horse heart protein, residue 39 for the yeast Candida krusei protein, and residue 62 for a triple mutant cytochrome c expressed via the yeast Saccharomyces cerevisiae. The structures of these protein variants were derived by Gray and co-workers by side-chain substitution of tuna and rice cytochrome c, whose X-ray crystal structures are kn0wn,~2a followed by molecular mechanics energy m i n i m i z a t i ~ n . ~The ~~,~ relative positions of the surface histidines are diagrammed in reference 2a. These structures were used without further modification to build the IAL Hamiltonian, which was then inverted to obtain the Green's function and thereby the superexchange couplings, or charge resonanceenergies, for the forward and back electron transfers in the three cytochrome c proteins,

given in Table VII. Provided for comparison are the values of the electronic factors extracted from the experimental rate data. The IAL method yields unprecedented agreement with experiment on the actual values of the superexchange charge resonance energies, A. The IAL values based on the forward 3ZnPorphyrin Ru(II1) value of the electronic energy E D A provide agreement to better than a factor of 2 with experiment. The sensitivity of the results to EDAis given two ways, as a finite difference derivative near EDA= -5.9 eV and near EDA= -6.8 eV and by comparison with the IAL results obtained at EDA= -6.8 eV for the back transfer. Note that the superexchangecharge resonance energies for the three proteins respond quitedifferently, with two of them changing sign; a behavior that highlights the interference characteristic of any full Hamiltonian nonperturbative treatment of superexchange. These IAL results do share the common property of any quantitative theoretical model that a shift of EDAdownward by 0.9 eV, all else being constant, will tend to enhance the rate of electron transfer mediated by hole superexchange, in particular by a factor of 13 in the case of His39 calculation.

-

-

IV. Discussion The matrix inversion approach provides a nonperturbative connection between the rate of quantum processes driven by indirect coupling and the corresponding electronic Hamiltonian and is drawn from the venerable Green's function formalism in scattering theory. Recent formal analyses in the electron-transfer context have provided algebraic elaborations based upon the perturbation expansion,43 upon iterative reduction of the propa g a t ~orr ~its~iterative b~ildup,2*9~~ and upon the reduction of the path summation to self-avoiding paths.'s Presented here is the direct noniterative numerical implementation in a realistic and practical form for molecular and macromolecular superexchange as demonstrated by the results presented in Tables V-VII. In particular the spectroscopic tests (PES data, Table V) do not depend on the accuracy of a prior factorization of the data into electronic and nuclear contributions and hence provides a clear test of the basic quantitative validity of this approach. For the optimum value of the through-bond parameter yo = 1.5 eV, all computed IAL values are within a factor of 3 of the PES splitting data and usually better than a factor of 2. This is better than one might anticipate given the schematic aperiodic lattice description of the bridge but provides evidence that the basic numerical parameters and the interaction energy function are clearly in the realistic range. Note that the interaction energy function connects all sites within 6 A, so that even these relatively covalent systems are not approximated as nearest-neighbor interaction through-bond chains. The IAL Hamiltonian provides distinctly superior results relative to the MNDO calculation for the orbital splitting, which undoubtedly reflects the general fact that a semiempirical method is best used for observables within

Model for Long-Range Electronic Superexchange the set (geometries,dipole moments, heats of formation) on which it was calibrated. The advantage of the IAL approach indeed emerges in the application to macromolecules, where the coarse-grained heavyatom-only site description is expected to be most appropriate, yet which enables the required matrix computations to proceed for a system of 5000 valence electrons in less than a minute on a vector supercomputer. The extended-Hiickel approach, which is superior for thecase of the intervalencetransitions of the binuclear dithiaspiro compounds, does not provide order of magnitude absolute accuracy when applied to proteins.33 It is likely that the realistic results obtained by the IAL Hamiltonian reflect the dominant through-space character of the three-dimensional network of couplings in the folded macromolecule. Extended HuckeP and IAL (Table V) can both perform well for dominantly covalent systems,though neither is as reliable (consistently better than a factor of 2) as would bedesirable for such smaller molecular systems. On the other hand, the IAL Hamiltonian contains an interaction energy function specifically calibrated against splitvalence ab initio results in the nonbonded range from van der Waals contact out to 6 A and performs far better in applications to proteins. Hence the IAL parameters, calibrated solely at the subunit (amino acid side-chain) level, evidently provide a more balanced treatment of both covalent and noncovalent contributions. The IAL strategy is the first method to provide factor of 2 or 3 accuracy in whole proteins. Among other advantages, the present analysis is then of sufficient realism to serve as a test and verification of our understanding of the mechanism of long-range electron transfer in proteins. Understanding the value of A at these distances of 13-16 A and in particular how it can be as large as experiment indicates, rather than as small as WKB would predict (HDAin Table VII; eq 2), is critical to quantitative insight into an element of protein design, namely, that proteins containing redox partners need not necessarily position them in direct proximity in order to obtain functional electron-transfer rates. The basic requirement of a successful and realistic theory of superexchange is that it quantitatively explains how the “insulating” medium intervening between the two special trap sites participates quantum mechanically, providing a coupling far worse than at contact but far better than separationthroughvacuum. The closed-shell electrons of the intervening amino acid residues are driven to delocalize by the off-diagonal exchange interactions but held tightly to their localized orbitals by the energy mismatch between their IPS and the EDA.The balance of these forces along interfering pathways throughout the entire macromolecular network is computed by the matrix inversion procedure, eq 3a, which then requiresaccurate values for the subunit diagonal energies (which governs the IPS) and for the off-diagonal elements (which govern the delocalization), and accurate values for the EDA. The emphasis here upon thorough calibration of the amino acid ionization energies and of the interaction energy function, and upon inclusion of all the valence orbitals of the cytochrome without prior perturbation pathway restriction, constitutes a proposal for the set of the most important characteristics in the theory of macromolecular electronic coupling. Many details at the subangstrom level are treated here only crudely, yet the remarkable agreement with experiment in Table VI1 is a significant indication that the indispensable characteristics are captured. The experimentally observed very weak but functional electronic delocalization is reproduced in the present quantum molecular model. The twotiered strategy is essential for accuracy, since the key subunit calibration step prevents approximations in the macromolecular representation of the valence electrons from deteriorating the accuracy of the important subunit energetic parameters. No special treatment has been given to hydrogen bonds, yet, even at this level of modeling, hydrogen bonds play a significant role in bringing atom-atom distances {N-O, N-N, 0-0) closer

The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 5591 than van der Waals contact, which then increasesthe off-diagonal element between their respective VBOs. However, due to the existence of the dense interpenetrating network of interfering pathways, a strengthening of a single off-diagonalmatrix element can decrease or increase the total charge resonance energy A, with nearly equal probability. Quantum interference effects must be present in any full Hamiltonian method, which renders the calculation of superexchange charge resonance energies more complex than single or few pathway perturbation calculations.53 Our whole protein matrix inversion description of electronic superexchange (eq 3) is equivalent to an (infinite) explicit path summation approach4* withsignedpath amplitudes. It is thenevident that theinterfering signs make it considerably more difficult to determine when a sum over a finite number of intertwined through-bond and through-space paths through a compact and folded macromolecule has converged. Furthermore analyses in progress of the participation of the various sites within the protein itldicate that all the orbitals in a broad zone (> 10 A cross section) between the redox sitesparticipate in determiningthe superexchangecoupling,which brings the utility of the intrinsically linear concept of pathways within a three-dimensional folded macromolecule into question. Due to the relative simplicity of a pathways approach, however, a careful examination of this is warranted, but it is to be emphasized that with the current full-protein numerical matrix approach in hand, any subdomain of proposed dominant importance can be studied systematically as an approximation to the whole matrix by extracting and inverting the correspondingmatrix sub-block. In this regard we have performed initial calculations on the basis of the Sherman-Morrison an important and exact sparse matrix method to rapidly determine the change in the inverse of a matrix, such as ( E i O - H)-l, given a change in the matrix, H.For example, the buried crystallographic water in the cytochrome c structure, H20 no. 26, was removed by this method, and the result in His39 protein electron transfer was to reduce A by 27%to +O. 145cm-I. Removal of all crystallographic waters yields +0.137 cm-I, hence the one water molecule is the notable participant. Removal of one &methyl group of Leu32, which is a critically placed residue for the His33 electron transfer, reduces the A markedly, to 4 . 0 6 7 cm-I for the CD1 methyl, whereas removal of the CD2 methyl is less significant, causing a 30% reduction in A to +0.13 cm-I. Subsequent protein conformational relaxation would lead to further changes, yet in the existing structure the CD1 methyl is vitally important, which implicates the critical participation of noncovalent interactions as hydrogenbonding is impossible. We are proceeding to use the ShermanMorrison formula as well as other analytical tools in further protein studies.

Conclusion In summary, we have described a minimal but quantitatively effective model Hamiltonian, which is used nonperturbatively to calculate the charge resonance energy upon which nonadiabatic long-range biomolecular electron-transfer rates depend. The successfulnumerical agreement obtained verifies the importance of detailed attention to the realistic quantitative representation of the inhomogeneous energetics of the protein medium, dictated by the sequence and ionization potentials, and to the parameters of the overall distance dependence used to assemble the offdiagonal elements of the Hamiltonian, obtained here by analysis of triple-c ab-initio results. The systematic calibration of the Hamiltonian is accomplished completely at the small molecule level with no adjustment to any macromolecular properties or data. Since no kinetic data have been used in the calibration, a fundamental conclusion is that hole superexchange as encapsulated in the IAL Hamiltonian is demonstrated to be kinetically competent as a mechanism to explain how charge resonance

5592

Gruschus and Kuki

The Journal of Physical Chemistry, Vol. 97, No. 21, 1993 Subunit Diagonal Energy Calibration

energies can be so significant despite large donor-acceptor separations. The IAL method can be used to estimate absolute values of the electronic charge resonance energies and not just their ratios. The inhomogeneous aperiodic lattice model is derived as a form of quantum molecular modeling appropriate for macromolecules rather than as a species of all-atom molecular orbital theory. As such it describes a relatively schematicrepresentation of the superexchangelattice, and considerable further improvement in predictive capability should be expected as more realistic side-chainwavefunctions,redox site orbitals, and detailed integrals are used to replace the elements of IAL. This systematic extension, based on valence bond orbitals, as well as further studies on cytochromec and its interactionswith other proteins, is in progress. Finally it should be emphasized that the computationally straightforward matrix inversion approach for the nonperturbative evaluation of electronic superexchange is directly generalizable to any Hamiltonian at any level of sophistication.

energy is calibrated from one member of the series. All data input into the calibration process are shown in square-cornered boxes in the flowchart. The valence bond orbitals (VBOs) are assigned by the rules shown in the flowchart to the heavy atoms by using the structure of the subunit analogue as a guide. In the case of proteins, the subunit analogue consists of a single a-amino acid whose backbone is “capped” by using backbone atoms of the residues before and after it to effectively serve as N-acetyl- and C-methylamidegroups. The alanine and glycine fragments of the easily- and non-easily-ionizable amino acids were first calibrated to the respective alanine and glycine IPS. In the second step the IAL Hamiltonian for the entire residue is constructed, with the orbital diagonal energy of the alanine or glycine component fixed at the predetermined value, and the diagonal energy parameter for sites in the functional group is then calibrated to reproduce upon diagonalization the experimental IP of the particular residue. For the non-easily-ionizable protein subunitcategory,the first IPScome from the amide orbital of the backbone rather than a side chain orbital as is the case with the easily ionizable residues. The atoms of the aliphatic side chains also number too few to use the obvious hydrocarbon analogues in the calibration process. In the PES of these residues, however, the onset of the u-bond molecular orbital energies, or the “u band”, is clear, so this energy was chosen to calibrate the orbital diagonal energies for the aliphatic side chains. The IAL Hamiltonian for the full medium consists of Hamiltonian subblocks for each subunit of the medium along with additional off-diagonal elements (eq 4) arising from intersubunit orbital interactions. The number of orbital sites for the cytochrome proteins is -2500. The inversion of (&A - H), the most time-consuming step in the computation of the superexchangecoupling A, takes about 30 s on an IBM ES/9000 supercomputer. In addition to cytochrome c, a 100 residue protein, we have performed IAL calculations on cytochrome c peroxidase, a -300 residue (-7000 valence orbital) protein. The maximum size system that can be calculated is limited by the matrix size allowed by the computer memory. Using the one gigabyte virtual memory of the IBM ES/9000,and assuming a symmetrical matrix, the maximum double-precision matrix is about 15000 X 15000, and for single precision 20000 X 20000. IAL Hamiltonians for proteins are sparse because of the 6-A cutoff for interactions; typically 90% or more of the matrix is zero, hence, sparsematrix routines could easily boost theallowable system size even further. The overall precision of the Lawdin partitioning method and the numerical IAL procedurewas tested by comparing analytically equivalent but numerically completely different formulas. In tests with N X N matrices from N = 110 to 2500, the numerical inversion yielded 13decimal place accuracy (out of 18calculated), indicating that the numerical routines used to calculate the IAL A are fully stable at this size.

Appendix

References and Notes

The IAL strategy is based upon dividing a macromolecule into subunitswhose properties are then calibrated from smallmolecule data or high-accuracy calculation. The orbital site diagonal energy of each new type of subunit (monomer unit, functional group, side chain, etc.) in the superexchange medium must be calibrated. Figure 5 shows the flow chart for the calibration process for easily- and non-easily-ionizable residues. The first step is determining the closest chemical analogue for which an experimental ionization potential (IP) or, in some cases, the photoelectron spectrum (PES) is available. For example, if adenosine were a part of the macromolecule, adenine could be used to calibrate the adenine subunit of the nucleoside. Subunit analogues with less than seven or eight heavy atoms should be avoided. In the special case of smaller molecules (Figures 2 and 3), the entire bridge becomes a single “subunit” and the diagonal

(1) Electron Transfer Reactions in Metalloproteins; Sigel, H., Sigel, A., Eds.;Marcel Dekker, Inc.: New York, 1991; Metal Ions in Biological Systems, Vol. 27. (2) (a) Wuttke,D. S.; Bjerrum, M. J.; Winkler, J. R.;Gray, H. B.Science 1992,256, 1007-1009. (b) Axup, A. W.; Albin, M.; Mayo, S. L.; Crutchley, R. J.; Gray, H. B. J . Am. Chem. SOC.1988, 110, 435-439. (c) Meade, T. J.; Gray, H. B.; Winkler, J. R. J . Am. Chem. SOC.1989, 1 1 1 , 4353-4356. (3) The Photosynthetic Bacterial Reaction Center: Structure and Dynamics; Breton, J., Vermeglio, A., Eds.; Plenum Press: New York, 1988. (4) Long-range Electron Transfer in Biology; Palmer, G.,Ed.; Springer-Verlag: Berlin, 1991; Structure and Bonding, Vol. 75. (5) Plato, M.; MBbius, K.; Michel-Beyerle, M. E.; Bixon, M.; Jortner, J. J . Am. Chem. SOC.1988,110, 1279-7285. (6) Marcus, R. A. Chem. Phys. Lett. 1988, 146, 13-22. (7) Miller, J. R.; Beitz, J. V.; Huddleston, R. K. J . Am. Chem. Soc. 1984, 106, 5057-5068. (8) Kuki, A.; Wolynes, P. G.Science 1987, 236, 1647-1652. (9) Kuki, A. In Long-Range Electron Transfer in Biology; Palmer, G. A., Ed.;Springer-Verlag: Berlin, 1991; Structure and Bonding, Vol. 75, pp 49-83.

I Protein method. (Small molecule bridge method is a simpliflcalion of the method beiow as the entire bridge is treated equivalently to a single residue.) For a given residue I, isolate the wordlnates plus the coordinates for the backbone C , 0.0nd Ca of i-1, and those for the backbone N and Ca of l i t

nuclear position (site). The number of orbitals is

of the side-chain sites 01 the residue. For rartdues

*

I

I

with easib ionizable side-chains. the alanine fracment in set to the diaoonal enerov cmde1ermIn.dlor lhe alanine fragment The c k t o m cdibrahon wlil deliver a diagonel energy for the v c m t side-chain unit. slerting at the gamme-atom. For aiiphatic residues, the glycine fragment is set to the diagonal energy predetermined lor glycine.

c

Hamiltonian for the given residue

1t,

\

(Construction of Interaction Function with Respect to Valence Bond Ortitai Overlap

I

\

Calculate the IAL interaction (Ea. 4) betwean all sites as a function of internuclear separation (zero

I

I

I

I 1

1 c

Diagonalize the-residue subblock Hemitonian

It calibraling Leu, 11% Val, Pro, *la, Arg, Lys residues with non.easilvionlzabie aide~chalns

c

I

c lonlxalion Polenllaia

P h o l o e i e s l r o n SpecIra

the IP, the calibration of the diagonal site energy is eigenvector with >10.h population (squared amplitude) in the non.

Figure 5. The flow chart for the calibration process for easily- and noneasily-ionizable residues.

-

Model for Long-Range Electronic Superexchange (10) Newton, M. D. Chem. Reu. 1991, 91, 767-792. (1 1) Schiff, L. I. Quantum Mechanics, 3rd ed.; McGraw-Hill, Inc.: New York, 1968; section 37. (12) Goldberger, M. L.; Watson, K. M. Collision Theory; Wiley: New York, 1964; section 6.5. (13) Weissbluth, M. Aroms and Molecules; Academic Press: New York, 1978; sections 9.4 and 14.4. (14) (a) Marcus, R. A. J. Chem. Phys. 1956,24,966. (b) Marcus, R. A. Annu. Rev. Phys. Chem. 1964, 15, 155. (15) Moser, C. C.; Keske, J. M.; Warncke, K.; Farid, R. S.; Dutton, P. L. Nature 1992, 355, 796-802. The average expotential falloff in proteins would correspond to a completely unrealistic 1.87-eV WKB barrier height if the simple one-particle tunneling formula (eq 2) were used. (16) Marcus, R. A.; Sutin, N. Biochim. Biophys. Acta 1985, 811, 265322. (17) Jortner, J. J . Chem. Phys. 1976, 64, 4860. (18) Hoptield, J. J. Proc. Narl. Acad. Sci. 1974, 71, 3640-3644. (19) Gamow, G.Z . Phys. 1928, 51, 204. (20) DeVault, D. Quantum Mechanical Tunnelingin BiologicalSystems, 2nd ed.; Cambridge University Press: Cambridge, 1984. (21) Gruschus, J.; Kuki, A. Chem. Phys. Lett. 1992, 192, 205-212. (22) Marchi, M.; Chandler, D. J . Chem. Phys. 1991, 95, 889-894. (23) Newman, W. H.; Kuki, A. J . Chem. Phys. 1992, 96, 1409-1417. (24) Hall, R. W.; Prince, M. R. J . Chem. Phys. 1991, 95, 5999-6004. (25) Onuchic, J. N.; Beratan, D. N. J . Chem. Phys. 1990,92,722-733. (26) Beratan, D. N.; Onuchic, J. N.; Betts, J. N.; Bowler, B. E.; Gray, H. B. J. Am. ChemSoc. 1990, 112, 7915-7921. (27) Therien, M. J.; Chang, J.; Raphael, A. L.; Bowler, B. E.; Gray, H. B. In Long-Range Electron Transfer in Biology; Palmer, G.,Ed.; Springer-Verlag: Berlin, 1991;Structure and Bonding, Vol. 75, pp 109-129. (28) Onuchic, J. N.; de Andrade, P. C. P.; Beratan, D. N. J. Chem. Phys. 1991, 95, 1131-1138. (29) Kuznetsov, A. M.; Ulstrup, J. J . Chem. Phys. 1981,75,2047-2054. (30) Ulstrup, J. Charge Transfer Processes in Condensed Media; Spring er: Berlin, 1979; Lecture Notes in Chemistry, Vol. 10. (31) Christensen, H. E. M.; Conrad, L. S.;Mikkelsen, K. V.; Nielsen, M. K.; Ulstrup, J. Inorg. Chem. 1990, 29, 2808-2816. (32) Broo, A,; Larsson, S.J . Phys. Chem. 1991, 95, 4925-4928. (33) Siddarth. P.: Marcus. R. A. J . Phvs. Chem. 1990. 94. 8430-8434. (34) Cave, R.'J.; Baxter, D. V.; W. A. Goddard, I.; Baldeschwieler,J. D. J. Chem. Phys. 1987,87,926-935. (35) Paddon-Row, M. N.; Engelhardt, L. M.; Skelton, B. W.; White, A. H.; Jsrgensen, F. S.;Patney, H. K. J. Chem. SOC.,Perkin Trans. II 1987, 1835. (36) Paddon-Row, M. N.; Jordan, K. D. In Molecular Strucrure and Energetics; Liebman, J. F., Greenberg, A., Eds.; VCH Publishers: 1990;Vol. 6, Chapter 6. (37) (a) Stein, C. A.; Taube, H. J . Am. Chem. SOC.1981,103,693-695. (b) Stein, C. A.; Lewis, N. A.; Seitz, G.J . Am. Chem. SOC.1982, 104,25962599. (38) Basu,G.; Kubasik, M.;Anglos, D.;Secor, B.; Kuki, A.J. Am. Chem. Soc. 1990, 112,9410-9411. (39) Closs, G.L.; Piotrowiak, J. M.; Fleming, G. R. J . Am. Chem. Soc. 1988, 110, 2652-2653. (40) Sigman, M. E.; Closs, G.L. J . Phys. Chem. 1991, 95, 5012-5017. (41) (a) SDarDaalione. M.: Mukamel, S.J . Phvs. Chem. 1987.91.3938. (b) Sparpaglibne; d.;Mukamel, S.J. Chem. Phis. 1988,88, 3263. (42) The standard perturbation pathway summation derives from the approximate expansion of the Lippmann-Schwinger equation, which is an alternative relation, similar to but distinct from eq 3a, involving the zerothorder Hamiltonian in the denominator? ) ) 43 To fulfill the present nonperturbative strategy, it is necessary to recast the implicit Lippmann-Schwinger equation as an explicit (matrix) equation, eq 3a, with the fully-coupled Hamiltonian in the denominator. (43) Ratner, M. A. J. Phys. Chem. 1990, 94, 4877-4883. (44) Evenson, J. W.; Karplus, M. J . Chem. Phys. 1992,96,5272-5278. (45) See also the analysis of the non-perturbative matrix inversion advantage in a different context, Section V of Maquet, A.; Chu, S.-I.; Reinhardt, W. P. Phys. Rev. A 1983, 27, 2946-2970. (46) The donor and acceptor trap sites are excluded from the H to be inverted. The purpose of the Green's function calculation is to capture the weak superexchange process at the moment the nuclear subsystem passes through the resonant transition state, and hence in the highly non-adiabatic limit the participation of the terminal donor and acceptor trap sites are held to minimum order in V. Further discussion of the indirect non-adiabatic rate theory will be presented elsewhere. In thecase of time-independentobservables

The Journal of Physical Chemistry, Vol. 97, No. 21. 1993 5593 such as PES, the condition that (E,O- f I - I q+,)has zero projection upon 14,) is equivalent to intermediate normalization.'" (47) Ostlund, A.; Szabo,N. S.Modern Quantum Chemistry. Introduction to Advanced Elecrronic Structure Theory; lst, revised ed.; McGraw-Hill: New York, 1989. (48) Sneddon, S.F.; Morgan, R. S.;Brooks, C. L. J . Phys. Chem. 1989, 93, 8115-8118. (49) Lijwdin, P.-0. J. Math. Phys. 1962, 3, 969. (50) Lbwdin, P.-0. J . Mol. Specrrosc. 1963, 10, 12. (51) Siddarth, P.; Marcus, R. A. J. Phys. Chem. 1990, 94, 2985-2989. (52) (a) Larsson,S. J. Am. ChemSoc. 1981,103,40344040. (b) Larsson, S.;Broo, A.; Kaellebring, B.; Volosov, A. Inr. J. Quantum Chem. Quantum Biol. Symp. 1988, 15, 1-22. (53) Liang, C.; Newton, M. D. J. Phys. Chem. 1992, 96, 2855-2866. (54) Baker, A. D.; Scharfman, R.; Stein, C. A. Tetrahedron Lett. 1983, 24,2957-2960. ( 5 5 ) The energies of unsolvated anion states of many relevant organic compounds are in the 0 to +2 eV range: Jordan, K. D.; Burrow, P. D. Acc. Chem. Res. 1978, 11, 341-348. Quasi-free. electrons in liquids have typical energies (VO)also near zero; *OS eV: Holroyd, R. A.; Tames, S.;Kennedy, A. J. Phys. Chem. 1975, 79, 2857-2861. (56) Newton's careful intercomparison of multielectron and independent electron results provides an important validation of the MO a p p r o a ~ h . ~ ~ . ~ ~ (57) Newton, M. D. J . Phys. Chem. 1988, 92, 3049-3056. (58) Yates, K. Hiickel Molecular Orbital Theory;Academic Press: New York, 1978. (59) Handy, N. C.; Marron, M. T.; Silverstone, H. J. Phys. Reu. 1969, 180,4548. (60) The values of 7 0 chosen here fall in the range expected for a throughbond coupling constant between adjacent valence bonds, as can be extracted by a valence bond analysis of Gaussian 6-jlG' molecular orbital energies of simple hydrocarbons. (61) Examination of various covalent hydrocarbon geometries indicated that pairs of the atom-centered VBOs used here have site separations d that overestimate the interaction distance between true bond-centered VBOs by roughly half a bond distance on average. A compensation of 0.75 A was hence substracted from all internuclear distances in the evaluation of V(d); this is achieved by setting a = tp + 0.75 A. (62) Alternatively, environmental effects upon the orbital energies can be included by including neighboring atoms in the subblock diagonalization; yet carrying out this more complicated calibration procedure gave the same IAL A values as the present method to within a factor of 2. (63) The non-easily-ionizable category also encompases trapped water molecules and the carbon atoms of the heme methyl, vinyl, and propionate side groups. (64) The backbone segments of all residues calibrated are 'capped" with an acetyl group and an N-methyl group taken from the backbone coordinates of the neighboring peptide bonds. (65) Cannington, P. H.; Ham, N. S. J. Electron Spectrosc. 1979, 15, 79-82. (66) Klasinc, L. J . Electron Spectrosc. 1976, 8, 161-164. (67) The amplitudes from orbital populations are 0.632 for Ru and 0.368 on the six NHj groups: Guenzburger, D.; Gamier, A.; Danon. J. Inorg. Chim. Acta 1977, 21, 119-131. (68) Vedeneyev, V. I.; Gurvich, L. V.; Kondrat'yev, V. N.; Medvedev, V. A.; Frankevich, Y. L. Bond Energies, Ionization Potentials, and Elecrron Affinities; Edward Arnold (Publishers) Ltd.: London, 1966. (69) Seddon, E. A.; Seddon, K. R. The Chemistry of Ruthenium; Elsevier: Amsterdam, 1984. (70) The El 2 values are measured in acetonitrile solution: Miller, L. L.; Nordblom, G. Mayeda, E. A. J. Org. Chem. 1972, 37, 916-918. (71) Stynes, H. C.; Ibers, J. A. Inorg. Chem. 1971, 10, 2304. (72) (a) Takano, T.; Dickerson, R. E. Proc. Nar. Acad. Sci. U.S.A.1980, 77,6371. (b) Therien, M. J.; Selman, M. S.; Chang, I.; Winkler, J. R.; Gray, H. B. J . Am. Chem. SOC.1990, 112, 2420-2422. (c) Bowler, B. E.; Meade, T. J.; Mayo, S.L.; Richards, J. H.; Gray, H. B. J . Am. Chem. SOC.1989,111, 8757-8759. (73) Magarshak, Y.; Malinsky, J.; Joran, A. D. J . Chem. Phys. 1991,95, 418-432. (74) da Gama, A. A. S. J. Theor. Biol. 1990, 142, 251-260. (75) Goldman, C. Phys. Rev. A 1991, 43, 4500-4509. (76) Ramsey, B. G. J . Org. Chem. 1979, 44, 2093-2097. (77) Matsumura-Inoue,T.; Kuroda, K.; Umezawa, Y.;Achiba, Y. J . Chem. SOC.,Faraday Trans. 2 1989,85, 857-866. (78) Press, W. H.; Flannery, B. P.; Teukolsky, S.A,; Vetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, 1986; p 66.

b.;