Competitive Electron Transfers in Model Triad Systems. MD Simulations

The dependence of the charge separation rate constant on the spatial arrangement of the subunits of the triad system, which arises from the present ca...
0 downloads 0 Views 860KB Size
Download PDF
J. Phys. Chem. 1996, 100, 8815-8826

8815

Competitive Electron Transfers in Model Triad Systems. MD Simulations Maria Hilczer*,†,‡ and M. Tachiya*,†,§ Department of Physical Chemistry, National Institute of Materials and Chemical Research, Tsukuba, Ibaraki 305, Japan, Department of Chemistry, The UniVersity of Tsukuba, Ibaraki 305, Japan, and Institute of Applied Radiation Chemistry, Technical UniVersity, Wro´ blewskiego 15, 93-590 Ło´ dz´, Poland ReceiVed: October 31, 1995; In Final Form: February 23, 1996X

The charge separation reactions for model systems containing an electron donor, D, and two electron acceptors, A, in acetonitrile solvent at room temperature have been investigated. The electron transfer rate constants for A1-D-A2 and D-A1-A2 systems with various angular arrangements of the subunits are expressed in terms of the two-dimensional statistical distributions of the reaction coordinates, which are defined as the electrostatic potential difference between the respective pair of the reacting sites. These distributions have been evaluated by the MD computer simulation method. The dependence of the charge separation rate constant on the spatial arrangement of the subunits of the triad system, which arises from the present calculations, has been compared with the corresponding result obtained on the basis of the dielectric continuum model of the solvent. The dynamics of the reaction coordinates for both model systems has been characterized by the classical time autocorrelation and cross-correlation functions.

I. Introduction Multichromophoric systems are recently the subject of intensive research, both experimentally1-25 and theoretically.26-39 Most of these studies are devoted to the investigation of photoinduced charge separation reactions as a means of capturing and storing solar energy. The capture of solar energy is extremely efficient in some biological systems, e.g., in the photosynthetic bacterial reaction centers. When the reaction center absorbs a photon, an electron is transferred from the excited bacteriochlorophyll dimer (BChl2*) to the about 17 Å separated bacteriopheophytin (BPh) in a very short time of approximately 3 ps. There is an additional bacteriochlorophyll molecule placed between BChl2 and BPh. Thus, this primary charge separation reaction actually occurs in the triad system. The role of the accessory BChl is, however, still unclear and the question of whether electron transfer occurs directly from BChl2 to BPh or through a distinct intermediate reduced state of BChl is extensively discussed in the literature.25-27 Understanding the photoinduced electron transfer reactions that occur in nature has an important practical implication. A long-term goal of this research is to design supramolecules that can mimic the charge separation function of the natural systems and, consequently, to control the conversion of solar energy and store it as useful chemical products. At present the number of studies devoted to the design and synthesis of complex supramolecular systems is continually increasing.1-24 Part of the synthesized supramolecules utilizes molecular compounds that are structurally reminiscent of those found in natural systems, e.g., porphyrins, quinones, carotenes, etc. Other systems are built from totally abiotic components like complexes of the second-row transition metals and small organic molecules. Several extensive review articles exist that cover the whole spectrum of synthesized supramolecular or multicomponent systems.1-10 * Correspondence can be sent to Professor Tachiya to National Institute of Materials and Chemical Research or to Dr. Hilczer to Technical University. † National Institute of Materials and Chemical Research. ‡ Present address: Technical University. § The University of Tsukuba. X Abstract published in AdVance ACS Abstracts, May 1, 1996.

S0022-3654(95)03213-8 CCC: $12.00

The synthetic models are mostly composed of three subunits,11-22 although some work has been done on tetrachromophoric23 and pentachromophoric24 systems. The artificial photosynthetic triads contain usually one donor and two acceptors (D-A1-A2 or A1-D-A2) or two donors plus one acceptor (D1-D2-A) subunits. The successive or competitive electron transfers that occur in supramolecular systems produce usually long-range chargeseparated states. The rate constants of these intramolecular processes are dependent on many factors, e.g., on donoracceptor distance and orientation, free energy change of reaction, electronic coupling between reacting centers, and also temperature and solvent. Van Dijk et al.20 studied recently some interesting trichromophoric systems in which the kinetics of electron transfer processes shows a very strong solvent dependence. The charge separation reactions in these systems can be controlled by small changes of solvent polarity, which has allowed the authors to produce one out of two or three different feasible charge transfer states. A large number of factors that influence the electron transfer processes in real systems make it impossible to control and investigate experimentally the importance of each of these factors separately. Many of the synthesized systems do not possess sufficient structural constraints to isolate the influence of such basic parameters as the donor-acceptor distance or orientation. The computer simulation techniques that are employed in studies of electron transfer reactions in supramolecules are free from these experimental limitations. In many cases they can provide a satisfactory explanation of the experimental observations on the molecular level. However, for biological systems, which are large and inherently complicated, the computer simulation is a time-consuming method and requires some simplifying approximations that can influence to a large extent the final obtained results. For instance, the shape and position of the free energy surfaces calculated by computer simulation for the primary electron transfers in a bacterial reaction center depend on the inclusion of electronic polarizability of the protein atoms and change considerably with the increase of the number of atoms incorporated explicitly into the computation. Particu© 1996 American Chemical Society

8816 J. Phys. Chem., Vol. 100, No. 21, 1996 larly, the role of solvent (water) molecules inside and around the protein-membrane system seems to be extremely important.34 In our opinion, simple model computer simulations of photoinitiated charge separation reactions in supramolecules can considerably deepen our understanding and explain how particular factors influence these processes. The present work has been stimulated by an increasing interest in the natural and synthesized many-chromophore supramolecules. A general goal of this paper is to present the correct approach to describe complex charge transfer processes in these systems. Let us consider the photosynthetic reaction center. There, electron transfers occur competitively between two paths, for example, from the special pair to the accessory chlorophyll or to the pheophytin, from the pheophytin to the quinone or back to the special pair, and so on. So far, these competitive electron transfers have been analyzed by using one-dimensional free energy curves. However, this is not appropriate, especially when the solvent dynamical effect is important. Here, we introduce two-dimensional free energy surfaces to treat competitive or consecutive electron transfers. For a simple model system we consider a triad system 1-2-3 embedded in a polar solvent and analyze the competition between electron transfers from 1 to 2 and from 1 to 3. We also discuss the results for another pair of charge separation reactions: the electron transfers from 2 to 1 and from 2 to 3. Our model triads are composed of one donor (D) and two acceptor (A1 and A2) centers. Thus, we employ the notation D-A1-A2 and A1-D-A2, respectively, to label the triad subunits 1, 2, and 3 that participate in the pairs of the electron transfer processes mentioned above. The triads are immersed in a polar, molecular solvent (acetonitrile) at room temperature. We are mostly interested in the influence of (a) solvent, (b) spatial arrangement of the triad subunits, and (c) the free energy changes on the competitive electron transfer processes that can occur in our model systems. We describe the electron transfer reactions between the photoexcited donor D* and one of the acceptors in terms of the reaction coordinates, e(VD - VA1) and e(VD - VA2), that represent the differences between the electrostatic potentials generated by solvent molecules at the donor and acceptor centers. As shown in refs 38 and 40, these coordinates are appropriate for studying electron transfers in the considered systems. It is important to note that the competitive (or successive) charge separation processes in a triad system are mutually dependent and haVe to be considered in the two-dimensional reaction coordinate space. The pair of coordinates, e(VD - VA1) and e(VD - VA2), is correlated, and as a measure of this correlation we can choose the correlation coefficient F or, for description of a dynamical behavior, the time cross-correlation functions. The coefficient F and the mean values mi and standard deviations σi characterize the two-dimensional statistical distribution φ of the reaction coordinates. The rate constants ki for electron transfers to acceptors A1 or A2 are proportional to the probability of finding the system on the intersection lines Γi of the free energy surfaces of the reactant and the respective product states. This probability is obtained by integrating distribution φ over Γi. Our definition of the reaction coordinates results in an extremely simple equation for Γi, i.e., the intersection line is a straight line parallel to one of the coordinate axes. Consequently, the calculation of the rate constants becomes easier than for other choices of the reaction coordinates.35 In the presentation that follows, we begin in section II with the description of details of the performed computer simulations.

Hilczer and Tachiya In section III we present the results for both systems, A1-DA2 and D-A1-A2, and discuss their dependence on the angular arrangement of the triad subunits. Particularly, we describe the stochastic properties of the reaction coordinates appropriate for competitive electron transfer processes, including the twodimensional and the marginal (one-dimensional) probability distribution functions and the dynamical correlations of the coordinates, which are expressed in terms of the time auto- and cross-correlation functions (section III.1). Next we calculate the rate constants for electron transfers to acceptors A1 or A2 together with the rate constant for the decay of the donor excited state and discuss their dependence on the free energy changes of both reactions (section III.2). These results are compared with the corresponding results obtained for the dielectric continuum model of the solvent. Finally, in section III.3, we discuss the structure of the solvent and the charge distribution in the vicinity of the triad, which result from the pair radial distribution functions of distances between the sites of acetonitrile molecules and the atoms of our supramolecular system. II. Method of Calculation The systems considered in the present calculations consist of three solute atoms labeled 1, 2, and 3 surrounded by 497 acetonitrile molecules. The distance R12 between the centers of atoms 1 and 2 is the same as the distance R23 between the centers of atoms 2 and 3 and equals 5.5 Å. The solutes are immobile and the distances R12 and R23 are kept constant during all simulation runs. The MD simulations are performed for several different angular arrangements of the solute atoms which are characterized by the angle β between the vectors R21 and R23 or, alternatively, by the angle γ between R12 and R13. The relation between these two angles for a given configuration of the solute atoms is γ ) (π - β)/2. The values of β assumed for particular simulation runs are 60, 90, 120, 150, and 180°. Two additional simulation runs were carried out, one for a pair of atoms separated by 5.5 Å and surrounded by 498 acetonitrile molecules and the other for the solute atom plus 499 solvent molecules. CH3CN-CH3CN interactions are described by the six-site potential developed by Bo¨hm et al.41 In this potential sites are located on the atoms of the acetonitrile molecule and each site of one molecule interacts with every site on another molecule through the standard 6-12 Lennard-Jones and Coulomb potentials. The acetonitrile molecule is treated as a rigid object, with bond lengths that are deduced from diffraction data on the liquid and equal to rC1-H ) 1.087 Å, rC1-C2 ) 1.460 Å, rC2-N ) 1.170 Å, rH-H ) 1.771 Å and the bond angle ∠(H-C1-C2) ) 109.8°. In the above notation C1 and C2 stand for the carbon atoms in the methyl and cyanide groups, respectively. This model potential is widely used in MD simulations and gives a good description of the structural, dynamical, and dielectric properties of both the bulk liquid42,43 and the acetonitrile clusters.44 The interactions of the solute atoms with the sites on a solvent molecule are restricted only to the 6-12 Lennard-Jones potential. A value of  ) 18kB is taken to be the Lennard-Jones well depth and s ) 4.9 Å to be the collision diameter for each of the solute atoms. The MD program used for the calculations is based on the MDMPOL coded by Smith and Fincham45 and calculates the time evolution of the system in the NVE ensemble using a leapfrog algorithm for the center-of-mass motions and a leapfrog-quaternion algorithm for the angular motions of the solvent molecules.46 The system of 500 particles was confined to a cubic box of side length 35.1214 Å, and the usual periodic boundary conditions and the Ewald summation of the electro-

Electron Transfer in Triad Systems static interactions were employed. The equations of motion were integrated with a time step of 2 fs, which ensured that the total energy was conserved within less than 0.08% over the course of a 80 ps trajectory in the equilibrium simulations. Equilibration of the system was performed for a period of 18 ps. During this period the linear and rotational kinetic energies of the system were rescaled to suit the assumed temperature of 291 K. In all equilibrium simulations the root-mean-square fluctuations of the resulting temperature were smaller than 4.9 K. A simulation run of about 450 ps was performed for each angular arrangement of the solute atoms. In order to avoid trapping of the solvent molecules into some artificial configurations in the vicinity of the solute atoms, this run was divided into five separate runs between which the system was first warmed up to a temperature of about 600 K and then gradually cooled back to 291 K. The electrostatic potential energies eV1, eV2, and eV3, produced by solvent molecules at each of the solute atoms, were calculated after every two time steps of the equilibrium trajectory of the system and were recorded sequentially during the production period of the simulation. Long-range electrostatic interactions were included in the potential experienced by the solute atoms by applying an anisotropic approximation to the Ewald summation.47 In this approximation the Ewald potential for a pair of point charges is given by the formula

1 E(r) ) + a1r2 + a2r4 + a3r6 + a4r8 + a5T4 + a6T6 + r a7T8 + a8T4r2 + a9T6r2 + a10T4r4 (1) where Tn ) xn + yn + zn and r ) (x2 + y2 + z2)1/2 denotes the distance between the test charge at the solute center and a partial charge on a particular atom of an acetonitrile molecule. Distances in eq 1 are expressed in units of side length of the simulation box, and the constants a1-a10 have the values 2.094 395, -4.506 792, 6.651 269, -10.866 13, 7.511 320, 17.071 59, 60.539 89, -23.279 44, -113.0078, and 65.196 80, respectively. Summation of the terms in eq 1, performed for all the partial charges of all the acetonitrile molecules in the simulation box, gives the electrostatic potential energy at the center of the solute atom. Formula 1 is relatively simple and allows us to obtain more precise values for the interaction energy than the spherical approximations to the Ewald sum.47-50 The electrostatic potential differences e∆V21 ) -e∆V12 ) e(V2 - V1), e∆V23 ) e(V2 - V3), and e∆V13 ) e(V1 - V3) between the respective pair of the solute atoms are random variables, which we denote as q1, q2, and q3, respectively. The random variables q1 and q2 can be used as the reaction coordinates for the competitive electron transfer processes in the model system A1-D-A2.38,40 Similarly, the random variables -q1 and q3 describe the electron transfers in the D-A1-A2 system. We constructed the histograms of the twodimensional probability distribution functions φ12(q1,q2) and φ13(-q1,q3) of the respective pairs of the reaction coordinates q1, q2 and -q1, q3 for all considered angular arrangements of the solute atoms using the data that were collected during the equilibrium simulation runs. Also, the marginal distribution functions φi(qi) of the particular random variables q1, q2, and q3 and the second-order characteristics, i.e., the correlation coefficients for the pairs of reaction coordinates, were calculated separately for each value of β. The performed simulation runs allowed us to obtain relatively smooth one-dimensional φi and two-dimensional φij probability distribution functions. They

J. Phys. Chem., Vol. 100, No. 21, 1996 8817 were also sufficient for calculation of the time autocorrelation and cross-correlation functions of the form

Ci(t) )

〈δqi(0) δqi(t)〉

Cc12(t) )

Cc13(t) )

〈δqi2〉

i ) 1, 2, 3

〈δq2(0)δq1(t)〉 F12(〈δq22〉〈δq12〉)1/2 〈δq3(0)δ[-q1(t)]〉 F13(〈δq32〉〈δq12〉)1/2

(2)

(3)

(4)

which characterize the classical dynamics of the considered random variables. In eqs 2-4 δqi(t) denotes the instantaneous fluctuation of qi(t) from its equilibrium average value, 〈δqi2〉 is the variance σi2 of the distribution φi(qi), and Fij is the coefficient of correlation between variables qi and qj. It is worth noting that two reaction coordinates qi and qj, which we choose to describe electron transfers in a given triad system, are directly connected to the third random variable qk. In the case of the model system A1-D-A2 the variable q3 ) q2 - q1, whereas for the D-A1-A2 triad we can write q2 ) q3 - (-q1). These equations together with basic formulas of probability theory provide some useful relations between the correlation coefficients Fij and the variances σi2 (or the standard deviations σi) of the reaction coordinates. For instance, the relation

σ32 ) σ12 + σ22 - 2σ1σ2F12

(5)

is valid in the system A1-D-A2, whereas for the D-A1-A2 model the analogous expression has the form

σ22 ) σ12 + σ32 - 2σ1σ3F13

(6)

Other simple relations

C3(t)σ32 ) C1(t)σ12 + C2(t)σ22 - 2σ1σ2Cc12(t)F12

(7)

C2(t)σ22 ) C1(t)σ12 + C3(t)σ32 - 2σ1σ3Cc13(t)F13

(8)

and

connect the time correlation functions Ci(t) and Cijc (t) for the reaction coordinates in the systems A1-D-A2 and D-A1A2, respectively. In the specific case of the triad composed of three identical atoms with the interatomic distances R12 ) R23, eqs 5-8 undergo considerable simplification and take the forms given in Appendix. The calculations were performed on an IBM RISC System/ 6000 Power Cluster, which executes 134 Mflops in a oneprocessor computation. III. Results and Discussion 1. Stochastic Properties of the Reaction Coordinates. A1-D-A2 System. Electron transfer can occur to one of the electron acceptor sites A1 or A2. We refer to the A1-D*-A2 state as redox state 0, A1--D+-A2 as state I, and A1-D+A2- as state II. The energy gaps ∆E0I and ∆E0II for the transitions 0 f I and 0 f II are expressed by the difference of the respective energies when the system is in the redox state 0, I, or II. Assuming that the electron interacts with its environment only electrostatically, these energy gaps can be expressed

8818 J. Phys. Chem., Vol. 100, No. 21, 1996

Hilczer and Tachiya

by the difference in the electrostatic energy of the transferred electron before and after transfer:

∆E0I ) e(V2 - V1) ) e∆V21

(9)

∆E0II ) e(V2 - V3) ) e∆V23

(10)

where ∆V21 ) V2 - V1 and ∆V23 ) V2 - V3 are the differences between electrostatic potentials generated at the donor D (atom 2) and the acceptor A1 (atom 1) and A2 (atom 3) positions, respectively. The energy differences

q1 ) e∆V21

(11)

q2 ) e∆V23

(12)

and

are, as shown before,38 the appropriate reaction coordinates for the description of the competitive electron transfer processes in the considered system. Within the framework of the linear response approximation both these coordinates are Gaussian random variables. Their joint probability distribution function φ12(q1,q2) is the well-known bivariate normal distribution

φ12(q1,q2) )

{

1

exp -

[( (

) ) ]}

q1 - m1 2 1 2 σ1 2(1 - F12)

2πσ1σ2x1 - F212 q1 - m1 q2 - m2 q2 - m2 2F12 + σ1 σ2 σ2

2

(13)

where m1 and m2 denote the equilibrium ensemble averages of the respective potential energy differences, e∆V21 and e∆V23, which are calculated for the redox state 0. Integrating φ12(q1,q2) over all values of q2 or q1 yields the marginal normal distributions φ1(q1) or φ2(q2) that are the probability functions for each of the reaction coordinates

φi(qi) )

1

x2πσi

{ ( )}

exp -

Figure 1. Distribution of the electrostatic potential energy difference between the donor and acceptor sites calculated for the symmetrical triad system A1-D-A2 in acetonitrile at 291 K. Results of the MD simulation (dotted line) are fitted by the Gaussian function (full line). Spatial arrangement of the triad subunits is described by the angle β ) 60° between the axes connecting D with A1 and D with A2.

1 qi - mi 2 σi

2

i ) 1, 2

(14)

The correlation coefficient F12, the standard deviations σ1 and σ2, and the mean values m1 and m2 of the random variables q1 and q2 are the only parameters of distributions 13 and 14. These parameters depend on the angular arrangement of the A1-DA2 system, and for each value of the angle β between axes D-A1 and D-A2 their values can be obtained from the electrostatic potential data collected during MD simulations. In the present calculations both acceptor sites are identical, and consequently, we have φ1(q1) ) φ2(q2) ≡ φ(q). Figure 1 shows an example of the marginal distribution of coordinate q for the A1-D-A2 system with β ) 60°. The dotted line corresponds to φ(q), which is evaluated directly from MD results, and the full line is the Gaussian approximation to φ(q) with the zero mean. Spatial arrangement of the A1-D-A2 system has a small influence on the standard deviation σ and on the mean value m of the electrostatic potential energy difference e∆V between donor and acceptor sites. For β ) 60° we obtained σ ) 0.335 eV, but for larger values of β the standard deviations were around 0.34 eV. In our calculation the parameters m(β) have small but nonzero values (on the order of 0.01 eV), which is caused by some differences in the solvation of the donor and acceptor sites. It will be explained in section 3 on the basis of the pair radial distribution functions of the distances between

Figure 2. Correlation coefficient F12 between reaction coordinates q1 and q2 as a function of the angle β in the A1-D-A2 triad system calculated from the MD simulation (black circles) and from the continuum model for two different values of the donor and acceptor radii: ra ) 2.45 Å (long dashed line) and ra ) 2.75 Å (solid line).

the subunits of our system and the particular atoms of acetonitrile molecules. It is worthwhile to note that the presence of the third atom does not change considerably the distribution of the potential difference. The parameters σ and m obtained from the simulation performed for a pair of solutes separated by R ) 5.5 Å in acetonitrile solvent are 0.34 and 0 eV, respectively. For comparison, the standard deviation for a pair of infinitely separated atoms is equal to 0.48 eV. This value was calculated to be (2)1/2σs, where σs is the standard deviation of the electrostatic potential at the center of a single solute atom surrounded by acetonitrile molecules. The dependence of the coefficient of correlation F12 between reaction coordinates q1 and q2 on the spatial configuration of the A1-D-A2 system is presented in Figure 2. In the figure we compare F12, which is obtained from MD simulations for various values of the angle β, with the corresponding results calculated on the basis of the dielectric continuum model. In the latter model the correlation coefficient is expressed by F12 ) λ12/(λ1λ2)1/2, where λ1 and λ2 denote the solvent reorganization energies associated with the reaction coordinates q1 and q2, respectively, and λ12 is an unnormalized measure of correlation between these coordinates.38 Data used for the continuum model calculations are as follows: temperature T ) 291 K, the optical dielectric constant op ) 1.81, and the static dielectric constant s ) 37.5. We present the dependence of σ12 on β for two different radii of the donor and acceptor sites: ra ) 2.45 Å and ra ) 2.75 Å. Both values are connected to the parameter σ of

Electron Transfer in Triad Systems

J. Phys. Chem., Vol. 100, No. 21, 1996 8819

Figure 3. Time correlation functions C(t) of fluctuations in the potential energy difference between donor and acceptor sites calculated for the symmetrical triad system A1-D-A2 in acetonitrile at 291 K for different values of the angle β.

solute atoms in the Lennard-Jones potential applied in the simulation. The former value is the hard core radius of the solute, ra ) s/2, and the latter corresponds to the minimum of the Lennard-Jones potential, ra ) 21/6s. The energies λ1, λ2, and λ12 were evaluated from the exact formula given in ref 51. Figure 2 shows that the coefficient of correlation between reaction coordinates q1 and q2 for a molecular model of solvent is nearly equal to F12 calculated on the basis of the continuum approximation for ra ) 2.75 Å for β in the range from 60° to about 100°. For larger values of β, including the linear configuration of the triad subunits, the correlation between reaction coordinates in the system A1-D-A2 is lower than that predicted from the continuum model. Figure 3 presents the equilibrium time correlation functions C(t) of fluctuations in the electrostatic potential difference between the donor and one of the acceptor sites (eq 2) calculated for various spatial arrangements of the A1-D-A2 system on the basis of MD simulations. All the correlation functions are very similar, especially for longer times (not shown in the figure). They are also similar to C(t) of the electrostatic potential difference obtained from the MD simulation of a pair of atoms49,50 as well as to the equilibrium solvation timecorrelation function calculated for a single solute atom in acetonitrile solvent.48,52 Like the last one, the correlation functions in Figure 3 have a bimodal character with a fast component very well described by a Gaussian and with a slow component that decays on a picosecond time scale. The molecular mechanisms responsible for this character of the solvation response are still a subject of discussion, but the fast part is usually interpreted as being due to small amplitude inertial motions of solvent molecules in the close vicinity of the solute. From the Gaussian fitting to the solvation response for very short times (e.g., from 0 to 0.05 ps), we can estimate the solvation frequency ω, related to t ) 0 curvature of C(t) by

ω2 )

( ) ∂2C(t) ∂t2

(15) t)0

The value of ω that results from the single atom simulation is equal to 11.6 ps-1. The characteristic frequencies ω, calculated for the acceptor atoms in the triad systems with different spatial arrangement β, are in the range 11.0-11.7 ps-1. For each value of β the value of ω for the donor site (central atom) is slightly smaller than for the acceptor atom and lies between 11.0 and 11.5 ps-1. The obtained frequencies imply a characteristic time for the fast solvation component of 85-91 fs. Similar values are reported for the fast component of the solvation response

Figure 4. Time cross-correlation functions Cc(t) of fluctuations in the potential energy difference between donor and acceptor sites calculated for the A1-D-A2 triad system for different values of the angle β.

for large solutes in acetonitrile solvent at room tempreature.52 The correlation functions C(t) in Figure 3 characterize the dynamics of the difference in the electrostatic potential between the donor and acceptor centers. Thus, the rate of decay of C(t) is related to the dynamics of acetonitrile molecules in the close vicinity of the triad. Small differences between C(t) observed for various spatial arrangements of the triad atoms can arise from different geometries of the space available for these solvent molecules. The characteristic frequency calculated from the time correlation functions presented in Figure 3 changes with β from 10.2 to 10.9 ps-1. Figure 4 shows the time cross-correlation functions (eq 3) for the reaction coordinates q1 and q2, calculated for various spatial arrangements of the triad subunits. The characteristic frequency of these functions decreases from ω ) 10.2 ps-1 for β ) 60° to ω ) 6.6 ps-1 for β ) 180°. This impies that the time constant for decay of the major Gaussian component of Cc(t) for the A1-D-A2 system increases from 98 to 152 fs when the arrangement of the triad atoms changes from a regular triangle to linear. D-A1-A2 System. In this system the D*-A1-A2 state is referred to as redox state 0 and D+-A1--A2 and D+-A1A2- are referred to as redox states I and II, respectively. The appropriate reaction coordinates for the description of the competitive electron transfer processes that can occur in the D-A1-A2 triad are

-q1 ) e(V1 - V2) ) eV12

(16)

q3 ) e(V1 - V3) ) eV13

(17)

where ∆V12 or ∆V13 are the differences between the electrostatic potential generated by acetonitrile molecules at the center of the donor (atom 1) and the acceptor A1 (atom 2) or the acceptor A2 (atom 3), respectively. Under the linear response approximation, both these coordinates are Gaussian random variables. Their joint probability distribution function φ13(-q1,q3) is the normal bivariate function

φ13(-q1,q3) ) 1 2πσ1σ3x1 - F132

{

exp -

[(

) ( ) ]}

q1 + m1 2 1 + σ1 2(1 - F132)

2F13

q1 + m1 q3 q3 + σ1 σ3 σ3

2

(18)

and their marginal distributions are univariate normal functions

8820 J. Phys. Chem., Vol. 100, No. 21, 1996

Hilczer and Tachiya

Figure 5. Distributions of the electrostatic potential energy difference between the donor and the acceptor A2 sites calculated for various spatial arrangements of the subunits of the triad system D-A1-A2 in acetonitrile at 291 K. The configuration of the triad is described by the angle γ between the axes D-A1 and D-A2.

of the form in (14) with i ) 1 or i ) 3. The mean of e∆V12 in the redox state 0 is -m1, whereas the average value of e∆V13 in this state equals zero for symmetry reason. The parameters of the Gaussian distributions of the reaction coordinates -q1 and q3 depend on the angular arrangement of the D-A1-A2 system that is described by the angle γ between the axes D-A1 and D-A2. Thus, the definition of the angle γ for the D-A1A2 system is analogous to the definition of the angle β in the case of the A1-D-A2 triad. The parameter σ1 of the Gaussian distribution of -q1 for the D-A1-A2 triad with a particular value of γ is, of course, equal to the standard deviation of the coordinate q1 for the A1-D-A2 system with the corresponding value of β. The Gaussian distributions of the reaction coordinate q3 ) e∆V13 for various values of γ are presented in Figure 5. The standard deviation σ3 changes from 0.432 to 0.335 eV when γ increases from 0 (linear arrangement of the triad) to 60°. The highest value of σ3 ) 0.437 eV was obtained at γ ) 15°. We can compare these values with σ3 calculated on the basis of the dielectric continuum model. It appears that, for each configuration of the triad, the distribution that arises from the latter model is narrower than the corresponding MD distribution. For a radius of the solute atoms ra ) 2.45 Å, the “continuum” value of σ3 decreases monotonously from 0.348 to 0.294 eV with an increase of the angle γ from 0 to 60°, and for ra ) 2.75 Å, the standard deviation changes from 0.323 to 0.264 eV over the same interval of γ. The variance of the distribution of the electrostatic potential difference between the donor D and the acceptor A2 sites, σ32(γ), is related to the reorganization energy λ3(γ) for the electron transfer process by

σ32(γ) ) 2kBTλ3(γ)

(19)

This equation allows us to calculate the dependence of the reorganization energy λ3 on the spatial configuration of the D-A1-A2 triad system. The reorganization energy λ3(γ) predicted from the MD simulations changes from 3.71 to 2.24 eV with an increase of the angle γ from 0 to 60°. As can be expected, the reorganization energy that results from the continuum model for both solute radii ra considered is smaller than the corresponding λ3(γ) predicted from the MD calculations. Figure 6 shows the dependence of the correletion coefficient F13 for the reaction coordinates -q1 and q3 on the angle γ between axes D-A1 and D-A2 in the triad D-A1-A2. For comparison we present also the functions F13(γ) calculated on

Figure 6. Correlation coefficient F13 between reaction coordinates -q1 and q3 as a function of the angle γ in the D-A1-A2 triad system calculated by the MD simulation (black circles) and by the continuum model for two different values of the donor and acceptor radii: ra ) 2.45 Å (long dashed line) and ra ) 2.75 Å (solid line).

Figure 7. Time correlation functions C(t) of fluctuations in the electrostatic potential energy difference between donor D and acceptor A2 sites calculated for the D-A1-A2 triad system in acetonitrile at 291 K for different values of the angle γ.

the basis of the continuum approximation for two different values of ra. As we can see, the correlation between coordinates -q1 and q3 obtained from the molecular model of the solvent is stronger than that predicted from the continuum model for both radii considered. The equilibrium time correlation functions of the reaction coordinate q3 calculated on the basis of MD simulations for various configurations of the D-A1-A2 system are plotted in Figure 7. The characteristic frequencies of the fast component of C3(t) are in the range 10.2-11.5 ps-1 for all considered arrangements of the triad subunits. At each γ the value of ω for the potential difference q3 is slightly larger than ω for the reaction coordinate q1 (cf. Figures 3 and 7). Thus, as we could expect, the time constant for decay of the fast Gaussian component of the response is longer for the potential difference between donor D and acceptor A1 than for the potential difference between D and A2. For each angular arrangement of the triad D-A1-A2 the time cross-correlation function Cc13(t) (eq 4) is shown to be equal to the corresponding C3(t) (see Appendix). By using the results in Figures 5 and 6, one can construct the joint probability distributions (eq 18) and therefore the twodimensional free energy surfaces that are presented in Figure 8. Figure 8a shows the free energy surface for the inital D*A1-A2 state in the case of the linear arrangement of the triad subunits. As we already pointed out, the electron transfer from D to A1 occurs when the potential difference -q1 takes a

Electron Transfer in Triad Systems

J. Phys. Chem., Vol. 100, No. 21, 1996 8821 triad) guarantees that the intersection ΓR for the pair of energy surfaces is described by a straight line that is parallel to one of the coordinate axes: qi ) qi# ) const and qj ∈ (-∞,∞).38 As mentioned before, such a simplicity is not ensured for other choices of the reaction coordinates for competitive electron transfer reactions in the triad systems.35 Within the framework of the linear response theory the joint probability density function, φij(qi,qj), is a two-dimensional Gaussian. Hence, the integral in eq 20 is reduced to the value φi(qi#) at qi ) qi# of the one-dimensional normal distribution of the reaction coordinate qi that is given by formula 14. A1-D-A2 System. In this system the intersection lines for the transfer processes 0 f I and 0 f II are q1 ) q1#, q2 ∈ (-∞,∞) and q1 ∈ (-∞,∞), q2 ) q2#, respectively. In the case of the charge separation reaction the values q1# and q2# for each spatial arrangement of the triad are given by

q1#(β) ) -λ1(β) - ∆G1

(21)

q2#(β) ) -λ2(β) - ∆G2

(22)

where ∆G1 and ∆G2 are the free energy changes of the corresponding reactions. The rate constants for the electron transfers from the donor D to the acceptors A1 or A2 are expressed by

k0fI(β) )

2π 2 J φ [q #(β)] p 1 1 1

(23)

k0fII(β) )

2π 2 J φ [q #(β)] p 2 2 2

(24)

and

Figure 8. Contour maps of the potential energy function for the locally excited D*-A1-A2 state of the triad system with the linear (a) and a regular triangle (b) arrangements of the subunits. The labels stand for the values of energy in units of electronvolts, and the difference between successive contours is 0.1 eV.

respectively. In eqs 23 and 24 Ji for i ) 1, 2 stands for the transfer integral53,54

specific value labelled as -q1#, which is determined by the electronic properties of the donor and the acceptor A1. On the other hand, the electron transfer from D to A2 occurs when q3 takes another specific value q3#. If the rate of the fluctuation of the surrounding solvent is sufficiently fast compared with the rate of electron transfer at -q1 ) -q1# and q3 ) q3#, the equilibrium distribution on this free energy surface is maintained. In this case the two electron transfers occur independently and each electron transfer can be treated using the respective one-dimensional free energy curve. However, if the rate of fluctuation of the solvent is not so fast, the two electron transfers interfere with each other. In this case one needs the two-dimensional free energy surface to analyze their competition. Figure 8b shows the free energy surface in the case of a regular triangle arrangement of the triad. As we can see, the shape of the free energy surface depends on the geometry of the supramolecule. 2. Electron Transfer Rates. The electron transfers between the redox states 0 f I and 0 f II of the triad system occur along the intersections Γ of the relevant free energy surfaces in the two-dimensional reaction coordinate space. The reaction rate constant kR for a given transition R is proportional to the probability of finding the system at ΓR

with R ) R12 ) 5.5 Å, and we take the constant parameters as a ) 1.0 Å-1, J0 ) 100 cm-1, and ra ) 2.75 Å. In the symmetrical A1-D-A2 triad the rate constant k0fI(β) is equal to k0fII(β), and because fluctuations of the solvent polarization are very fast, the rates of these two processes are additive. Thus, the rate constant for the decay of the donor excited state due to both electron transfer reactions is equal to k(β) ) k0fI(β) + k0fII(β). In the MD calculations the reorganization energy λi depends slightly on the arrangement of the triad atoms, and consequently, the value of k changes with β. The variation of k is, however, rather small. For example, for ∆G1 ) ∆G2 ) -1.0 eV, the rate constant is equal to 3.7 ns-1 at β ) 60° and decreases to 2.0 ns-1 when the triad configuration becomes linear. The reorganization energy λi expressed by the approximate equation of the dielectric continuum model does not depend on the angular configuration of the triad. Consequently, the rate k is only a function of the size of the solutes, ra, and the free energies of both reactions, ∆Gi for i ) 1, 2. Taking ra ) 2.75 Å and ∆G1 ) ∆G2 ) -0.8 eV, we obtain k ) 0.4 ps-1. For lower ∆G ) -1.0 eV, we have k ) 1.6 ps-1. It is worth noting that when we calculate the reorganization energy from the continuum model in a more sophisticated way, i.e., by taking into account the two- and three-center correction terms,51 we obtain values of λ and k that show small angular variation. However, since an increase of k0 f I is observed only in the very close vicinity of β ) 60° for the contact pair separations of the triad subunits, we will hereafter compare the

kR ≈ ∫∫ φij(qi,qj) dqidqj ΓR

(20)

where R ) 0 f I or 0 f II. Our choice of the reaction coordinates (q1, q2 for A1-D-A2 and -q1, q3 for D-A1-A2

Ji ) J0 exp[-R(R - 2ra)/2]

(25)

8822 J. Phys. Chem., Vol. 100, No. 21, 1996

Hilczer and Tachiya

Figure 9. Electron transfer rate constants k0 f I (lines 1 and 1′) and k0 f II (lines 2 and 2′) in the D-A1-A2 system in acetonitrile solution as functions of the angular arrangement of the triad subunits obtained for the continuum model (curves 1′ and 2′) and for the molecular model of the solvent (curves 1 and 2). Curves 3 and 3′ show the rate constants for the decay of the donor excited state. Calculations were performed for ∆G1 ) -0.8 eV, ∆G3 ) -1.6 eV, and ra ) 2.75 Å.

Figure 10. Electron transfer rate constant k0 f II in the D-A1-A2 system in acetonitrile solution as a function of the angular arrangement of the triad subunits, calculated based on the molecular model of the solvent for a variety of values of the free energy change ∆G3 of the reaction 0 f II. The label of each curve stands for the value of ∆G3 in units of electronvolts.

MD rate constants with the rate constants that arise from the approximate formulas of the continuum model. D-A1-A2 System. Here, the intersection line of the energy surfaces 0 and I is

-q1 ) -q1#(γ) ) -λ1(γ) - ∆G1 and the rate constant k0

f I(γ)

k0fI(γ) )

(26)

equals

2π 2 J φ [-q1#(γ)] p 1 1

(27)

The intersection line for 0 f II transfer is

q3 ) q3#(γ) ) -λ3(γ) - ∆G3 and the rate constant k0

f II(γ)

k0fII(γ) )

(28)

is given by

2π 2 J φ [q #(γ)] p 3 3 3

(29)

with

J3(γ) ) J0 exp[-a(2R cos γ - 2ra)/2]

(30)

In Figure 9 we present the rate constants for two competitive electron transfer reactions and the rate k(γ) ) k0 f I(γ) + k0 f II(γ) for the decay of the excited state of the donor as functions of the angular arrangement of the D-A1-A2 triad for both the molecular and the continuum model of the solvent. The free energy changes for electron transfers to the acceptors A1 and A2 were assumed to be -0.8 and -1.6 eV, respectively. As we can see, the γ-dependence of all the rate constants is stronger in the case of MD results. For instance, k0 f II increases from 28 ms-1 to 0.29 ps-1 when the configuration of the system changes from linear to a regular triangle. The same change of the angle γ causes the increase of the 0 f II rate constant only by a factor of 103 for the continuum approximation. In the range of γ from 0 to about 40° the decay of the donor excited state D* is governed exclusively by the electron transfer to acceptor A1. The 0 f II transition becomes important for γ larger than 45°. Figure 10 shows the angular dependence of the rate constant k0 f II for a variety of values of the free energy change, ∆G3, calculated for the D-A1-A2 triad in acetonitrile solution on the basis of MD simulation. The reorganization energy λ3(γ) connected with the 0 f II transition is over the whole range of

Figure 11. Electron transfer rate constant k0 f II in the D-A1-A2 system in acetonitrile solution as a function of the angular arrangement of the triad subunits, calculated based on the continuum model of the solvent for a variety of values of the free energy change ∆G3 of the reaction 0 f II. The label of each curve stands for the value of ∆G3 in units of electronvolts, and ra ) 2.75 Å.

γ larger than -∆G3 only for the cases of the free energy change equal to -1.6 and -2.0 eV. For lower values of ∆G3 the rates k0 f II(γ) partly enter the Marcus inverted region. A shallow minimum that appears in each of the curves results from the fact that the standard deviation and, consequently, the reorganization enegry λ obtained by MD simulations are slightly higher for γ ) 15° than for any other considered configuration of the triad.55 Figure 11 presents the k0 f II(γ) functions calculated for various values of the free energy change for the dielectric continuum approximation. Since the reorganization energy of the continuum model is lower than λ3(γ) for the molecular model of the solvent, the curve with ∆G3 ) -1.6 eV is partly, and the curves with lower ∆G3 are entirely, in the inverted region. The k0 f II decreases quite rapidly with an increase of γ for low values of the free energy change, e.g., the rate constant that is calculated for ∆G3 ) -3.2 eV equals 14.8 µs-1 for a linear system and decreases to 77.9 s-1 for a regular triangle arrangement of the triad subunits. The difference between the corresponding pair of curves that are plotted in Figures 10 and 11 for the same value of the free energy change results partly from different relative magnitudes of ∆G3 and λ3 for each curve of the pair. The curve k0 f II(γ) that arises from the molecular model for a given value of ∆G3 resembles the curve of the continuum model for higher free

Electron Transfer in Triad Systems

Figure 12. Rate constant for the decay of the excited donor state in the system D*-A1-A2 in acetonitrile as a function of the angular arrangement of the triad subunits, calculated based on the MD results for a variety of values of the free energy change. The label of each curve stands for the value of ∆G1 in units of electronvolts. For each curve ∆G3 ) 2∆G1.

energy change. This can be seen from a comparison of the curve labeled -1.6 in Figure 10 and the one labeled -0.7 in Figure 11. The difference between the reorganization energy λ3(60°) and ∆G3 for the continuum model with ∆G3 ) -0.7 eV is equal to that for the molecular model with ∆G3 ) -1.6 eV. It is necessary to note, however, that the variation of λ3 over the whole range of γ in the MD calculations is more than twice as large as that in the continuum approach. Thus, even if we shift the scale of ∆G3 as described above, the corresponding pair of curves considered above still gives ratios k0 f II(60°)/k0 f II(0°) that are different by a factor of 116. On the other hand, in the inverted Marcus region the rate k0 f II(γ) that is calculated on the basis of MD simulations decreases over the whole range of γ faster than k0 f II(γ) of the continuum approach. We can obtain a more quantitative relation by comparison of the MD rate constant for ∆G3 ) -4.8 eV (not shown in Figure 10) with the curve labeled -3.2 in Figure 11. For these values of the free energy change the respective sums λ3(0°) + ∆G3 for both models are equal to -1.1 eV. In this case the ratio k0 f II(60°)/k0 f II(0°) is 1.18 × 10-9 for the molecular model and 5.26 × 10-6 for the continuum approximation. Figure 12 shows the decay rate k of the excited donor state as a function of the spatial configuration of the triad system D*-A1-A2 calculated on the basis of MD simulations for a variety of values of the free energy change. For higher values of ∆G3 the curves are similar to curve 3 in Figure 9. When the value of ∆G3 gradually decreases, the γ-dependence of the rate constant k is more and more flat, which means that the transition 0 f I becomes dominating over the transition 0 f II over the whole range of γ. For extremely low values of the free energy change (not shown in the figure) the donor excited state decays at the rate described exclusively by k0 f I(γ) independently of the arrangement of the triad subunits. The values of the free energy change in Figures 9-12 are evidently too low in comparison with ∆G for any real triad system. For instance, the experimental value of ∆G3 in the photosynthetic reaction center is on the order of -0.26 eV.35 The system considered in the present paper is, however, a model system that is composed of three atoms in acetonitrile solution. Our simulation is not devoted to considering the real structure of solutes. Moreover, our calculations do not take into account the effect of the electronic polarization of the solvent molecules. This effect can change considerably the fluctuations in the electrostatic potential difference and, consequently, the calculated reorganization energy and ∆G connected with electron

J. Phys. Chem., Vol. 100, No. 21, 1996 8823 transfer reactions.36 The influence of the electronic polarization was introduced a posteriori in some MD simulations of primary electron transfer in the photosynthetic reaction center32,35 by dividing the calculated reorganization energy by the optical dielectric constant of the medium op ≈ 1.9. This approach is, however, criticized by Warshel et al.36 as being incorrect in the case of the charge separation reaction. It is claimed by other workers that inclusion of the high-frequency modes causes an increase of the reorganization energy (at least for nonadiabatic electron transfer reactions in highly polar solvents)56 and reduction of the transition rate.57 3. Solvent Structure around the Triad System. Within the framework of the dielectric continuum model the mean values, m1 ) 〈e∆V21〉 and m2 ) 〈e∆V23〉, of the electrostatic potential difference between the central atom 2 and one of the side atoms (1 or 3) are both equal to zero regardless of the angular arrangement of the triad. If we employ, however, molecular models of the solvent surrounding the triad, this property of the reaction coordinates distribution is not generally valid. The model applied in the present paper results in positive values of the mean electrostatic potential energies 〈eVi〉 produced by acetonitrile molecules at all the triad atoms. Moreover, the value of 〈eV2〉 for each configuration (except a regular triangle) of the triad is slightly higher than the average potential energies 〈eV1〉 and 〈eV3〉. Consequently, the mean values of the potential energy differences, m1 ) 〈eV2〉 - 〈eV1〉 ) 〈eV2〉 - 〈eV3〉 ) m2, are positive. The highest m1 ) 0.01 eV was obtained for the configuration with β ) 120° and the lowest value of 0.003 eV for the linear triad. Since these values are very small in comparison with the average energies 〈eVi〉, we do not attempt to give any quantitative angular dependence. However, it is worthwhile to note that the solute-solvent radial distribution functions (RDFs) calculated from our model for each atom of the triad and each site of the acetonitrile molecule show a clear difference in solvation of the central atom and the terminal atoms. In Figure 13 we present a set of the atom-site RDFs, g(r), calculated for the spatial configuration of the triad defined by the angle β ) 120° (or, equivalently, by γ ) 30°). The functions for the C1, C2, and N sites of acetonitrile molecules have only one distinct peaksthis shows that the correlation of solvent molecules around each solute atom is limited to the first solvation shell. These correlations for all the triad atoms are weaker than in the case of a single atom immersed in acetonitrile. For example, the ratios of the peak heights of RDFs for atoms 1 or 3 in Figure 13 (curves labeled 1) to those of RDFs for a single atom are equal to 0.88, 0.93, and 0.93 for C1, C2, and N sites, respectively. There is no distinct peak in the atom-hydrogen RDFs for the triad. Figure 13 also shows that the heights of the peak of the distributions g(r) for the central atom (curves labeled 2) are lower than those of RDFs for subunits 1 or 3. The functions g(r) alone are insufficient to determine the average configuration of acetonitrile molecules in the vicinity of each of the triad atoms. We can point out, however, some features characteristic for the first solvation layer of the triad. Since the site C2 is close to the center of mass of CH3CN, we first consider the function g(r) for the atom-C2 distance distribution. The minimum of this RDF for the side atoms occurs at r ) 6.65 Å, and the integrated coordination number up to this point is nC2 ) 11.2. The minimum of the central atom RDF is at the same value of r, but it gives nC2 ) 9.9. Moreover, the mean numbers of C1, N, and H sites in a sphere of radius r ) 6.65 Å around the terminal atoms are 11.0, 10.95, and 32.8, respectively, while the corresponding numbers for the central atom are 10.0, 9.7, and 30.4. Thus, the side atoms of

8824 J. Phys. Chem., Vol. 100, No. 21, 1996

Hilczer and Tachiya

(a)

(b)

(c)

(d)

Figure 13. Radial distribution functions describing the distributions of distances between each of the side atoms (curve 1) or the central atom (curve 2) of the triad and the particular sites of the acetonitrile molecules. C1 and C2 stand for carbon sites of the methyl and cyanide groups, respectively.

the triad compared to the central atom are solvated, on the average, by one more acetonitrile molecule. The number of solvent molecules in the first solvation shells of the triad atoms is approximately equal to the coordination number for bulk acetonitrile that is defined in terms of center-of-mass separations of the molecules.41 However, the spatial extent of the first solvation layer is shorter in pure solvent (the first minimum on the C2-C2 RDF occurs at r ) 6.15 Å). The main peak in the C2-terminal atom g(r) is at 5.05 Å, and it is shifted by 0.65 and 0.8 Å toward larger r in comparison with the main peaks in the RDFs for N and C1 sites, respectively. It means that the acetonitrile molecules in the close vicinity of the terminal atoms prefer the orientations in which their nitrogens or CH3 groups are directed toward the atom. The detailed analysis of the radial distribution functions in Figure 13 and the integrated coordination numbers allowed us to predict the average orientations of approximately eight acetonitrile molecules around each of the terminal triad atoms. We have found that five molecules are oriented with the nitrogen directed toward the atom (four with the average angle ∠(C2N-atom) about 90° and one with 170°) and the other three with the CH3 group toward the atom (two with ∠(C2-C1-atom) close to 100° and one with 130°). From the results shown in Figure 13 we cannot say, however, how these molecules are mutually oriented. We can only conjecture that the former molecules have their nitrogens directed toward the protons of the latter molecules. A similar analysis performed for the central atom of the triad shows that there are five solvent molecules with the nitrogen directed toward the atom, two with the average ∠(C2-N-atom) of about 96°, two with 130°, and one with 170°. Two or three other solvent molecules with CH3 toward the atom have the average angle ∠(C2-C1-atom) of 110° and 140°.

Figure 14. Total charge of the sites of the acetonitrile molecules, which are enclosed within a sphere of radius r around each of the side atoms (curve 1) and the central atom (curve 2) of the triad, as a function of r. Curve 3 shows the difference between curves 2 and 1.

We can calculate the total charge T(r) inside a sphere of radius r around each of the triad centers by evaluating the integrated coordination numbers nk(r) at the distance r for all sites of the acetonitrile molecule and summing nk(r) multiplied by the respective partial charges. The dependence of the total charge on the radius r that we have obtained for the terminal and the central atoms of the triad for β ) 120° is shown in Figure 14. The closest to the solute atoms are hydrogens of CH3CN. The coordination number for H sites reaches nH ) 1 at a distance r ) 3.4 Å from the side atoms or at r ) 3.5 Å from the central atom. At these distances the values of nk calculated from respective RDFs for other sites of CH3CN are very small, e.g., nN is equal to 0.0067 and 0.024 for the terminal and the central atoms, respectively. Since the partial charge of the H site is positive, the total charge T(r) in the close vicinity of all the

Electron Transfer in Triad Systems triad centers is also positive. Curve 3 in the figure shows the difference of the total charge functions T(r) for the central and the side triad units. In the region r < 4.5 Å, T(r) is larger for the terminal atoms by less than 0.07 e, but in the next region 4.5 < r < 7 Å, the difference in T(r) is positive with values up to 0.1 e. These mutual relations of the charge distributions around the triad atoms result in the slightly positive average values, m21 and m23, of the difference in electrostatic potential energy between atom 2 and one of the terminal atoms. The RDFs for other arrangements of the solute system are similar to the functions presented in Figure 13. For all considered values of the angle β the corresponding peaks of RDFs for distances between the terminal units and the C1, N, and H sites of acetonitrile molecules have the same position and intensity. When the angle β increases, the first peak of the side atom-C2 RDF becomes split into two separate maxima. In the case of the linear system these maxima occur at 4.55 and 5.15 Å and are of nearly the same height. The number of acetonitrile molecules in the first solvation shell of the triad atoms does not clearly depend on β. The only exception is the central atom of the triad with a rectangular configuration for which the coordination number, defined by the position of the first minimum of the C2-atom RDF (r ) 6.3 Å), is smaller than for other β and equals 7.9. The functions T(r) of Figure 14 are also representative of other angular arrangements of the traid centers. The subtle interplay of the positive and negative values of these functions for the central and the side units causes a weak angular dependence of the mean values of the reaction coordinates q1 and q2. IV. Concluding Remarks We have considered the reactions of electron transfer from photoexcited donor D* to one of the acceptor centers, A1 or A2, in the model triad systems, A1-D*-A2 and D*-A1A2, immersed in acetonitrile solvent at room temperature. To treat these competitive electron transfers, we have introduced the two-dimensional free energy surface Π0(qi,qj) for a given triad whose shape depends on the angular arrangement of the reacting centers (cf. Figure 8). Each reaction occurs at the intersection of the free energy surfaces of the reactant and the respective product, which in our approach is the straight line parallel to one of the coordinate axes qj or qi. The electron transfer rate constants have been expressed in terms of the twodimensional statistical distributions φij(qi,qj) of both solvent polarization coordinates. These distributions for the traid systems of various geometries have been evaluated by MD computer simulations. The statistical description of the reaction coordinates and the angular dependence of the electron transfer rate constants, which arise from the molecular model of solvent, differ considerably from predictions of the dielectric continuum model. The latter model underestimates the fluctuations of the electrostatic potential in the triad vicinity. Consequently, the statistical distribution φij(qi,qj) and the marginal distributions φi(qi) of the reaction coordinates calculated on the basis of MD results are wider than the respective “continuum” distributions. Moreover, the molecular model results in a stronger dependence of the correlation coefficient between coordinates qi and qj on the geometry of the supramolecule. As a simple model system for the photosynthetic reaction center, several groups are trying to construct triad systems in which back electron transfer is reduced. In order to reduce back electron transfer, what electronic properties should the components of the triad have? What geometry should the supramolecule have, and what static and dynamic properties should the

J. Phys. Chem., Vol. 100, No. 21, 1996 8825

Figure 15. Contour map of the potential energy function for the ionic triad system D+-A1--A2 with the linear arrangement of the subunits. The forward electron transfer occurs along the line q2 ) q2# and the backward electron transfer along the line q1 ) q1#. The labels stand for the values of energy in units of electronvolts, and the difference between successive contours is 0.05 eV.

solvent have? As illustrated in Figure 15, forward and back electron transfers in the D+-A1--A2 triad system occur along the lines q2 ) q2# and q1 ) q1#, respectively. Changing the electronic properties of the triad components corresponds to changing the positions of these lines. Changing the geometry of the supramolecule or the static property of the solvent corresponds to changing the shape of the free energy surface. Changing the dynamical property of the solvent corresponds to changing the rate at which the system fluctuates on the surface. On the basis of this kind of figure, we can optimize the electronic properties of the components, the geometry of the supramolecule, and the static and dynamical properties of the solvent to reduce back electron transfer. Thus, we can give useful guidelines for designing triad systems. Appendix If the triad is composed of three identical subunits and both terminal units are equally distant from the central unit, then the relations connecting the statistical parameters of the reaction coordinates q1 (or -q1), q2, and q3 have the forms

|F13| )

σ3 2σ1

(A1)

F12 ) 1 - 2F132

(A2)

2F132C3(t) ) C1(t) - F12Cc12(t)

(A3)

C3(t) ) Cc13(t)

(A4)

Equations A1 and A2 are fulfilled by the ensemble-averaged values of the standard deviations σi and the correlation coefficients Fij calculated for our triad systems, A1-D-A2 or D-A1-A2, for any angular arrangement of the atoms. Equations A3 and A4 connect the time autocorrelation functions Ci(t) (eq 2) and the time cross-correlation functions Ccij(t) (eqs 3 and 4). The last expression concerns the pair of the reaction coordinates, -q1 and q3, that is related to the competitive electron transfer processes in the D-A1-A2 triad. It shows that autocorrelation of the electrostatic potential difference between the donor D and the acceptor A2 decays at the same

8826 J. Phys. Chem., Vol. 100, No. 21, 1996

Hilczer and Tachiya

rate as cross-correlation of the potential differences between the donor and the acceptor A2 and between the donor and the acceptor A1. Taking into account eqs 2, 4, A1, and A4, we obtain

2〈δq3(0)δ[-q1(t)]〉 ) 〈δq3(0)δq3(t)〉

(A5)

which in the notation of the probability theory has the form

2 cov[q3(0),-q1(t)] ) var[q3(0),q3(t)]

(A6)

where cov[q3(0),-q1(t)] is the covariance function of the random variables q3 and -q1 and var[q3(0),q3(t)] is the variance function of the random variable q3. Thus, the correlation between q3(0) and the value of the reaction coordinate -q1 after a time interval t is, on the average, half as strong as the correlation between values q3(0) and q3(t). Relation A6 is valid for all times t, and at t ) 0 it is equivalent to eq A1. References and Notes (1) Wasielewski, M. R. Chem. ReV. 1992, 92, 435. (b) Wasielewski, M. R. In Photoinduced Electron Transfer; Fox, M. A., Chanon, M., Eds.; Elsevier: Amsterdam, 1988; Part A, p 161. (2) Gust, D.; Moore, T. A.; Moore, A. L. In AlternatiVe Fuels and the EnVironment; Sterrett, F. S., Ed.; Lewis: Boca Raton, FL, 1995; p 125. (3) Gust, D.; Moore, T. A. Top. Curr. Chem. 1991, 159, 103. (4) Gust, D.; Moore, T. A. AdV. Photochem. 1991, 16, 1. (b) Gust, D.; Moore, T. A. Science 1989, 244, 35. (c) Gust, D.; Moore, T. A. In Supramolecular Photochemistry; Balzani, V., Ed.; Reidel: Dordrecht, 1987; p 267. (5) Connolly, J. S.; Bolton, J. R. In Photoinduced Electron Transfer; Fox, M. A., Chanon, M., Eds.; Elsevier: Amsterdam, 1988; Part D, p 303. (6) Balzani, V. Tetrahedron 1992, 48, 10443. (7) Sauvage, J.-P.; Collin, J.-P.; Chambron, J.-C.; Guillerez, S.; Coudret, C.; Balzani, V.; Barigelletti, F.; De Cola, L.; Flamigni, L. Chem. ReV. 1994, 94, 993. (8) Balzani, V.; Campagna, S.; Denti, G.; Serroni, S. In Photoprocesses in Transition Metal Complexes. Biosystems and Other Molecules. Experiment and Theory; Kochanski, E., Ed.; Kluver: Dordrecht, 1992; p 233. (9) Denti, G.; Campagna, S.; Sabatino, L.; Serroni, S.; Ciano, M.; Balzani, V. In Photochemical ConVersion and Storage of Solar Energy; Pelizzetti, E., Schiavello, M., Eds.; Kluwer: Dordrecht, 1991; p 27. (10) Sabbatini, N.; Guardigli, M.; Lehn, J. M. Coord. Chem. ReV. 1993, 123, 201. (11) Hung, S.-C.; Macpherson, A. N.; Lin, S.; Liddell, P. A.; Seely, G. R.; Moore, A. L.; Moore, T. A.; Gust, D. J. Am. Chem. Soc. 1995, 117, 1657. Hung, S.-C.; Lin, S.; Macpherson, A. N.; DeGraziano, J. M.; Kerrigan, P. K.; Liddell, P. A.; Moore, A. L.; Moore, T. A.; Gust, D J. Photochem. Photobiol. A 1994, 77, 207. (12) Gust, D.; Moore, T. A.; Moore, A. L.; Gao, F.; Luttrull, D.; DeGraziano, J. M.; Ma, X. C.; Makings, L. R.; Lee, S.-J.; Trier, T. T.; Bittersmann, E.; Seely, G. R.; Woodward, S.; Bensasson, R. V.; Rougee, M.; De Schryver, F. C.; Van der Auweraer, M. J. Am. Chem. Soc. 1991, 113, 3638. (13) Johnson, D. G.; Niemczyk, M. P.; Minsek, D. W.; Wiederrecht, G. P.; Svec, W. A.; Gaines, G. L.; Wasielewski, M. R. J. Am. Chem. Soc. 1993, 115, 5692. (14) Wiederrecht, G. P.; Watanabe, S.; Wasielewski, M. R. Chem. Phys. 1993, 176, 601. (15) Harriman, A.; Odobel, F.; Sauvage, J.-P. J. Am. Chem. Soc. 1994, 116, 5481. (16) Collin, J.-P.; Harriman, A.; Heitz, V.; Odobel, F.; Sauvage, J.-P. J. Am. Chem. Soc. 1994, 116, 5679. (17) Brouwer, A. M.; Mout, R. D.; Maassen van den Brink, P. H.; van Ramesdonk, H. J.; Verhoeven, J. W.; Jonker, S. A.; Warman, J. M. Chem. Phys. Lett. 1991, 186, 481. Brouwer, A. M.; Mout, R. D.; Maassen van den Brink, P. H.; van Ramesdonk, H. J.; Verhoeven, J. W.; Warman, J. M.; Jonker, S. A. Chem. Phys. Lett. 1991, 180, 556. (18) Brouwer, A. M.; Eijckelhoff, C.; Willemse, R. J., Verhoeven, J. W.; Schuddeboom, W.; Warman, J. M. J. Am. Chem. Soc. 1993, 115, 2988. (19) Lawson, J. M.; Paddon-Row, M. N.; Schuddeboom, W.; Warman, J. M.; Clayton, A. H. A.; Ghiggino, K. P. J. Phys. Chem. 1993, 97, 13099. (20) van Dijk, S. I.; Wiering, P. G.; Groen, C. P.; Brouwer, A. M.; Verhoeven, J. W.; Schuddeboom, W.; Warman, J. M. J. Chem. Soc., Faraday Trans. 1995, 91, 2107. (21) Osuka, A.; Nagata, T.; Kobayashi, F.; Zhang, R. P.; Maruyama, K.; Mataga, N.; Asahi, T.; Ohno, T.; Nozaki, K. Chem. Phys. Lett. 1992, 199, 302.

(22) Osuka, A.; Marumo, S.; Maruyama, K.; Mataga, N.; Tanaka, Y.; Taniguchi, S.; Okada, T.; Yamazaki, I.; Nishimura, Y. Bull. Chem. Soc. Jpn. 1995, 68, 262. (23) Osuka, A.; Nakajima, S.; Maruyama, K.; Mataga, N.; Asahi, T.; Yamazaki, I.; Nishimura, Y.; Ohno, T.; Nozaki, K. J. Am. Chem. Soc. 1993, 115, 4577. Lee, S.-J.; DeGraziano, J. M.; Macpherson, A. N.; Shin, E.-J.; Kerrigan, P. K.; Seely, G. R.; Moore, A. L.; Moore, T. A.; Gust, D. Chem. Phys. 1993, 176, 321. (24) Gust, D.; Moore, T. A.; Moore, A. L.; Macpherson, A. N.; Lopez, A.; DeGraziano, J. M.; Gouni, I.; Bittersmann, E.; Seely, G. R.; Gao, F.; Nieman, R. A.; Ma, X. C.; Demanche, L. J.; Hung, S.-Ch.; Luttrull, D. K.; Lee, S.-J.; Kerrigan, P. K. J. Am. Chem. Soc. 1993, 115, 11141. Gust, D.; Moore, T. A.; Moore, A. L.; Lee, S.-J.; Bittersmann, E.; Luttrull, D. K.; Rehms, A. A.; DeGraziano, J. M.; Ma, X. C.; Gao, F.; Belford, R. E.; Trier, T. T. Science 1990, 248, 199. (25) Small, G. J. Chem. Phys. 1995, 197, 239. (26) Skourtis, S. S.; Mukamel, S. Chem. Phys. 1995, 197, 367. (27) Bixon, M.; Jortner, J.; Michel-Beyerle, M. E. Chem. Phys. 1995, 197, 389. (28) Thompson, M. A.; Schenter, G. K. J. Phys. Chem. 1995, 99, 6374. (29) Thompson, M. A.; Zerner, M. C. J. Am. Chem. Soc. 1991, 113, 8210. (30) Warshel, A.; Parson, W. W. Annu. ReV. Phys. Chem. 1991, 42, 279. (31) Bixon, M.; Jortner, J. J. Phys Chem. 1986, 90, 3795. (32) Schulten, K.; Tesch, M. Chem. Phys. 1991, 158, 421. (33) Creighton, S.; Hwang, J.-K.; Warshel, A.; Parson, W. W.; Norris, J. Biochemistry 1988, 27, 774. (34) Parson, W. W.; Chu, Z.-T.; Warshel, A. Biochim. Biophys. Acta 1990, 1017, 251. (35) Marchi, M.; Gehlen, J. N.; Chandler, D.; Newton, M. J. Am. Chem. Soc. 1993, 115, 4178. (36) Warshel, A.; Chu, Z. T.; Parson, W. W. J. Photochem. Photobiol. A 1994, 82, 123. (37) Chandler, D.; Gehlen, J. N.; Marchi, M. In Ultrafast Reaction Dynamics and SolVent Effects; Gauduel, Y., Rossky, P. J., Eds.; American Institute of Physics: New York, 1994; p 50. (38) Najbar, J.; Tachiya, M. J. Phys. Chem. 1994, 98, 199. (39) Fushiki, M.; Tachiya, M. J. Phys. Chem. 1994, 98, 10762. (40) Tachiya, M. J. Phys. Chem. 1993, 97, 5911. (41) Bo¨hm, H. J.; McDonald, I. R.; Madden, P. A. Mol. Phys. 1983, 49, 347. (42) Bo¨hm, H. J.; Lynden-Bell, R. M.; Madden, P. A.; McDonald, I. R. Mol. Phys. 1984, 51, 761. (43) Ohba, T.; Ikawa, S. Mol. Phys. 1991, 73, 999. (44) Al-Mubarak, A. S.; Del Mistro, G.; Lethbridge, P. G.; Abdul-Sattar, N. Y.; Stace, A. J. Faraday Discuss. Chem. Soc. 1988, 86, 209. (45) Smith, W.; Fincham, D. CCP5 Program Library; Science and Research Council, Daresbury Laboratory: Warrington, UK, 1982. (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987; p 196. (47) Adams, D. J.; Dubey, G. S. J. Comput. Phys. 1987, 72, 156. (48) Maroncelli, M. J. Chem. Phys. 1991, 94, 2084. (49) Hilczer, M.; Tachiya, M. J. Mol. Liq. 1995, 64, 113. (50) Tachiya, M.; Hilczer, M. In Ultrafast Reaction Dynamics and SolVent Effects; Gauduel, Y., Rossky, P. J., Eds.; American Institute of Physics: New York, 1994; p 447. (51) Motylewski, T.; Najbar, J.; Tachiya, M. To be published. (52) Maroncelli, M.; Kumar, P. V.; Papazyan, A.; Horng, M. L.; Rosenthal, S. J.; Fleming, G. R. In Ultrafast Reaction Dynamics and SolVent Effects; Gauduel, Y., Rossky, P. J., Eds.; American Institute of Physics: New York, 1994; p 310. (53) Newton, M. D.; Sutin, N. Annu. ReV. Phys. Chem. 1984, 35, 437. (54) Miller, J. R. In Supramolecular Photochemistry; Balzani, V., Ed.; Reidel: Dordrecht, 1987. (55) Calculating k0 f II(γ) for the D-A1-A2 triad, we have approximated the values of the reorganization energy λ obtained by MD simulations for several values of γ by a polynomial of the form λ(γ) ) 3.736 75 + 0.008 16γ - 0.000 55γ2. The nonlinear regression has been performed using the standard Levenberg-Marquardt algorithm with the goodness of fit defined by the coefficient of determination r2 ) 0.9976. The flat maximum of this polynomial within the considered interval of γ is, however, not at γ ) 15° but at γ ) 7.4°. (56) Kim, H. J.; Hynes, J. T. J. Phys. Chem. 1990, 94, 2736. (57) Zhu, J.; Cukier, R. I. J. Chem Phys. 1995, 102, 8398.

JP953213X