J . Phys. Chem. 1984, 88, 6152-6157
6152
(Bayer)), b (Ti02 (Bayer)-2.7% Li), and c (Ti02 (Bayer)-7% Li), the pHo values found were 6.5, 7.5 and 9.0, respectively. From the pHo values obtained in eq 7, and with the knowledge that the MV2+/MV+ redox potential is -0.69 V vs. SCE, Ef values of -0.23, -0.19, and -0.13 V are obtained for the compounds shown in traces a-c. Li doping shifts pHo in eq 7 to more basic pH values.
2ol
Figure 9. Dependence of rate of the change of photocurrent with time (AilAt) on pH for stirred and N2-purgedTiOzsuspension photocell. Cell conditions: TiOzpowder (250 mg); H,O (100 mL); [NaOAc] = 1.0 M; [KN03] = 0.1 M; [MVz+]= 1 mM. Platinum collector electrode at -0.20 V vs. SCE: (a) Ti02 (Bayer), (b) TiO, (Bayer)-2.7% Li, (c) TiO,
(Bayer)-7% Li. the Li-doped semiconductor powders used. The source of this effect is the shift of the Fermi energy level Ef with pH as shown by eq 7. Extrapolation of the two lines Ai/At when large and Ef = Ef(pH 0) - 0.059pH (at 25 "C) (7) small Ai/At changes were observed afforded a precision of fO.l pH unit in the determination of pH,. The intersection point (pH,) represents the pH value in eq 7 when the Fermi level (Ef) equals the redox potential of the couple MV2+/MVC. For traces a (Ti02
Conclusion In conclusion, an active Li-doped anatase catalyst in water photocleavage has been presented in this study. The first evidence has been presented here in which Li doping modifies the dynamics of charge transfer across the anatase interface, improving its water photocleavage efficiency. The link between the nature of the catalysis and the role of Li ions has been presented by using diverse physicochemical techniques. We think that the effect of doping on the electronic transfer in semiconductor powders can, to a considerable degree, be explained by the phenomena described in this paper. The role of surface hydroxyl groups in H, generation and O2photoadsorption is of great importance in photocatalytic reactions as confirmed by pH dependency of the products observed. This observation has been previously reported by Somorjai.1',31 Evidence has been presented which shows that modification of the semiconductor bulk by Li doping also has a profound effect on the photocatalytic products observed. The detailed mechanism of water photocleavage remains to be identified.
Acknowledgment. This study was supported by the Swiss National Science Foundation. We thank Dr. Panek of Bayer AG for the preparation of the catalyst samples and M. Buffat for his assistance with the electron microscopy work. Helpful discussions with Professor M. Gratzel are appreciated. C.M. thanks the European Photochemical Association for a travel grant. Registry No. H2, 1333-74-0;Li, 7439-93-2;water, 7732-18-5;anatase, 1317-70-0. (31) Ferrer, S.; Somorjai, G. J. Phys. Chem. 1981, 85, 1464.
Diffusion-Controlled Reactions between a Spherical Target and Dumbell Dimer by Brownian Dynamics Slmulation S. A. Allison,*+ N. Srinivasan, J. A. McCammon,* Department of Chemistry, University of Houston, University Park, Houston, Texas 77004
and S . H. Northrup Department of Chemistry, Tennessee Technological University, Cookeville. Tennessee 38505 (Received: June 18, 1984)
A new Brownian dynamics trajectory approach used recently to study the diffusion-controlled reaction of spherical reactants is extended to the simplest case of structured reactants: dumbell dimers reacting with a spherical target. It is shown that, for dimers with a single reactive subunit, electrostatic torques exerted on the dimer by the target can increase the reaction rate by "steering" the dimer toward productive collision geometries. The effects of variations in the reactive surface of the dimer and in the Coulombic and hydrodynamic interactions between the reactants are also considered.
1. Introduction since the pioneering work of Smo~uchowskiand Debye, who investigated the problem of diffusion-controlled reactions in the absence' and presence2 of centrosymmetric Coulombic forces, there has been a proliferation of theoretical studies of such reactions based on more complex models. These have included the effects +Current address: Department of Chemistry, Georgia State University, University Plaza, Atlanta, GA 30303.
0022-3654/84/2088-6152$01.50/0
of hydrodynamic interaction's4 and solvent caging effects.5 The general area of enzyme-catalyzed bimolecular reactions has received special attention. The rate of diffusional encounter of an enzyme and ligand may be influenced by many factors, such as (1) Smoluchowski, M. V. Phys. Z . 1916, 17, 557. (2) Debye, P. Trans. Electrochem. SOC.1942, 82, 265. (3) Friedman, H. L. J . Phys. Chem. 1966, 70, 3931. (4) Deutch, J. M.; Felderhof, B. U. J . Chem. Phys. 1973, 59, 1669. (5) Northrup, S.H.; Hynes, J. T. J . Chem. Phys. 1979, 71, 871.
0 1984 American Chemical Society
Diffusion-Controlled Target-Dimer Reactions
The Journal of Physical Chemistry. Vol. 88, No. 25, 1984 6153
orientation dependence of reactivity for one or both specie^,"^ internal-configuration-dependent reactivity,’@’2interactions with the protein surface outside the active ~ i t e , ’ ~noncentrosymmetric ?’~ electrostatic forces,I5 and enzyme saturation,16 to name a few. It has been proposed that recognition of the proper substrate by enzymes actually begins at the diffusional encounter stage.” These and other factors are reviewed in greater detail e l ~ e w h e r e . ’ ~ J ~ - ’ ~ When a variety of interactions are operating simultaneously, e.g., when detailed models are used for an enzyme and substrate, there is, unfortunately, little hope of obtaining an analytical solution for the rate constant. Numerical treatments then become necessary. One numerical method that has been used for solving the diffusion equation with orientational anisotropy effects combines finite differences in the radial coordinate and eigenfunction expansions in orientation space.20 A more recent technique, which is used in this work, is the Brownian dynamics trajectory method.21 In Brownian dynamics, diffusional trajectories of the system are computed stochastically by a series of small displacements. The Figure 1. Model system. The target sphere (subunit 1) is uniformly rate constant and mechanistic information are obtained from an reactive with radius a, = 2 A. The dimer consists of two touching spheres analysis of the trajectories. This method is sufficiently general (subunits 2 and 3) with radii a2 = a3 = 0.5 A. In most but not all to model systems of arbitrary configurational complexity, arbitrary simulations, only subunit 2 of the dimer is reactive. Charges q,, q2, and inter- and intramolecular forces, and allows for inclusion of hyq3 are placed at the centers of the respective subunits. Also shown are drodynamic interactions. the relative coordinates R,r, and 0 as well as b, the starting distance for the simulation and q, which defines the truncation sphere. The purpose of the present work is to study the combined effects of hydrodynamic interaction, anisotropic reactivity, and interThe Si vector components represent stochastic displacements due particle forces on diffusion-controlled bimolecular reactions for to collisions with the solvent and are obtained with the multivariate a simple model system in which one has structured reactants. In Gaussian random number generator GGNSM from the IMSL the next section, we outline the Brownian dynamics trajectory subroutine library.23 Di;”is the initial hydrodynamic interaction method and present the model system used in this work. In section tensor between subunits i and j . In this work, hydrodynamic 3, the method of relating collision probabilities to bimolecular rate interaction is approximated by the Oseen tensor with slip boundary constants is reviewed and extended. In sections 4 and 5 the details conditions. This representation has been shown to provide a of the simulation and the results are discussed. The two quantities reasonable and simple point force description of the relative upon which we concentrate are the actual bimolecular rate condiffusion of finite spheres at small separation^.^^ In this case, stants and a conditional probability related to the orientation of the diffusion tensor is given by the reacting species when they are near each other. We conclude with a summary of the results.
2. Brownian Dynamics To simulate the Brownian motion of a system of interacting particles modeled as beads or arrays of beads, we used the method of Ermak and McCammon.22 If the initial position of subunit i is r; in a space-fixed reference frame, its position after an unconstrained dynamics time step of duration At is
where
T j is the Oseen tensor
(S,Sj) = 2Di;”At (6) Solc, K.; Stockmayer, W. H. Int. J . Chem. Kinet. 1973, 5, 733. (7) Schurr, J. M.; Schmitz, K. S. J . Phys. Chem. 1976, 80, 1934. (8) Samson, R.; Deutch, J. M. J . Chem. Phys. 1978, 68, 285. (9) Shoup, D.; Lipari, G.; Szabo, A. Eiophys. J . 1981, 36, 697. (10) McCammon, J. A.; Northrup, S. H. Nature (London)1981,293, 316. (1 1) Szabo, A.; Shoup, D.; Northrup, S. H.; McCammon, J. A. J. Chem. Phys. 1982, 77, 4484. (12) McCammon, J. A. Rep. Progr. Phys. 1984, 47, 1. (13) Chou, K. C.; Zhou, G. P. J . Am. Chem. SOC.1982, 104, 1409. (14) Zhou, G.; Wong, M. T.; Zhou, G. Q.Eiophys. Chem. 1983,18, 125. (15) van Leeuwen, J. W. FEBS Lert. 1983, 156, 262. (16) Berg, 0. G.; Ehrenberg, M. Eiophys. Chem. 1983, 17, 13. (17) Neumann, E. In “Structural and Functional Aspects of Enzyme Catalysis”; Eggerer, H., Huber,R., Eds.; Springer-Verlag: West Berlin, 198 1; pp 45-48. (18) Noyes, R. M. Prog. React. Kinet. 1961, 1, 128. (19) Calef, D. F.; Deutch, J. M. Annu. Rev. Phys. Chem. 1983,34,493. (20) Zientra, G. P.; Nagy, J. A.; Freed, J. H. J . Chem. Phys. 1980, 73, 5092. (21) Northrup, S. H.; Allison, S. A,; McCammon, J. A. J . Chem. Phys. 1984,80, 1517. (22) Ermak, D. L.; McCammon, J. A. J . Chem. Phys. 1978, 69, 1352.
(3)
(4) a = ai
where kBis the Boltzmann constant, Tis the absolute temperature, and F;” is the initial force acting on subunit j excluding stochastic (solvent) forces and forces of constraint. Si is a vector of Gaussian random numbers of zero mean and variance-covariance
1
I + (1 - hij)Tj
+ aj
- rij
rij C ai rij 2 ai
+ aj
+ a,
(5)
and a, is the radius of subunit i . The model used is shown in Figure 1. The target sphere (subunit 1) is taken to be uniformly reactive with a radius of 2 A. The dimer consists of two touching spheres (subunits 2 and 3) each with a radius of 0.5 A. In most simulations, only one of the dimer spheres (subunit 2) is reactive. Forces between the target and dimer arise from charges q l , q2,and q3 (in electronic charge units, e = 4.8 X esu) which lie at the centers of spheres 1, 2, and 3, respectively. The important relative coordinates are R = 1/2(r2+ r3) - rl, r = r3 - r2, and B, the angle between R and r. From eq 1 and the above definitions of the relative coordinates, we have
At r’ = ro + -C(D3;” - D2jO).Fj0 kBT i
+ (S3- S2)
(7)
(23) “IMSL Library 7 Reference Manual”; IMSL International and Statistical Libraries, Inc.: Houston, TX, 1975. (24) Wolynes, P. G.; Deutch, J. M. J . Chem. Phys. 1976, 65, 450, 2030.
6154
The Journal of Physical Chemistry, Vol. 88, No. 25, 1984
The reason for the primes in eq 6 and 7 is that forces of constraint have not yet been accounted for. In the present problem, the distance between the two subunits of the dimer is fixed. As in previous work, subunit displacement is a two-step process.25 First, an unconstrained dynamics step is carried out (eq 6 and 7). Second, displacement correction vectors are added to correct for the constraint (see Appendix of ref 25). The corrected displacements turn out to be
R = R'
+I 200( Y ) ( D l 2 O - Dl3O).ro r = r'
(&ar^)
+ y2
- ro
At -(Ddo kB
+ D1l)*Fdo
&(At)
(9)
(10)
+ +
where Ddo = '/2(D22 D230)is the dimer translational diffusion tensor, Fdo= F20 F30is the net direct force acting on the dimer, and ( SRSR) = 2(Dd0 D I ,)At. The corresponding equation for r can be written
+
r = ro + 0 x ro At 0 = -D,(F20 2k~T
- F30) X ro + S,
to bimolecular rate constants is presented. As depicted in Figure 1, a sphere of radius b divides the relative separation space, R, into an outer region (R > b ) and an inner region (R < b). The value of b is chosen sufficiently large so that (1) interparticle forces (direct and hydrodynamic) are approximately centrosymmetric for R > b, and (2) the ensemble reactive flux through the R = b surface should have no angle dependence. Under steady-state conditions, the actual diffusion-controlled rate constant is given by2'
(8)
where d is the fixed distance between the dimer subunits and Do = kBT/(4rqa2)is the (slip) translational diffusion constant of a sphere of radius u2. It should be emphasized that eq 6-9 are only valid for short time steps. The time step, Ar, must be sufficiently short so that direct forces remain essentially constant during Ar and that 1 8- r^I b. Since there remains a finite probability that a pair achieving separation R = q would ultimately react, it is necessary to correct the measured recombination probability, /3, to account for this. For the special case where all reactive surface collisions lead to a reaction
= kD(b)/kD(q)
(15)
The general case where all reactive surface collisions do not necessarily lead to reaction is treated in ref 21. Equation 14 relates /3, obtained from the simulation, to the bimolecular rate constant, k. That part of the problem which involves diffusion in the vicinity of a target of arbitrary complexity is determined numerically (p). The remainder of the problem, which involves relative diffusion at large distances where the interactions are effectively centrosymmetric, is determined analytically (kD(b)). For particles interacting via a centrosymmetric potential of mean force, U(r), the following result may be employedS
(12)
where X denotes the vector cross product, D, = 3kBT/(32~qa23) is the slip rotational diffusion constant of the dimer for reorientation of its long axis, and (S,S,) = 2ID,At. Equation 10 corresponds to a pure translational displacement and eq 11 and 12 to a pure rotation of r. The first term on the right-hand side of eq 12 represents rotation of the dimer due to direct forces and the second term represents rotation due to Brownian motion.
3. Calculation of Rate Constants In the previous section, a method was outlined for computing how a system evolves in time during a Brownian dynamics simulation. This method has been used successfully to study internal motions in polymer^,^^^^^ transport coefficients of both rigid and semiflexible assemblies of spheres,25 and diffusion-influenced biomolecular reactions.21 For diffusion-influenced biomolecular reactions, one is ordinarily most interested in obtaining a bimolecular rate constant in order to make with experimental studies. To obtain a rate constant by Brownian dynamics one would, in principle, need to simulate a large ensemble of reactant pairs diffusing from large separations to the reaction surface. A more convenient quantity to measure by Brownian dynamics simulation is the recombination probability /3 for a pair of reactants diffusing in a finite domain. In ref 21, a formal derivation connecting /3 (25) Allison, S. A.; McCammon, .I. A. Biopolymers 1984, 23, 167. (26) Pear, M. R.;McCammon, J. A. J . Chem. Phys. 1981, 74, 6922. (27) Allison, S. A,; McCammon, J. A. Biopolymers 1984, 23, 363.
where D(R) is the relative diffusion constant. This result is not strictly valid for most cases studied here since the potential of mean force for motion of the dimer relative to the target sphere has an angular (e) dependence. Referring to Figure 1, let Qe be the charge on subunit 1 and AQ'e be the charge on subunits 3 and 2, respectively. Boltzmann averaging over 0 leads (by a lengthy derivation) to an equation of the same form as (16) with U(R)
N
--( P)~( ~ B TdQ'Qe2 6 ckBTR2
1
-)
+ -1 Dd
20 DM
(17)
where d = 2a2 is the distance between the centers of the two dimer subunits, E is the dielectric constant of the medium, and
D(R)=DM+-
Dd
120
(
dQ'Qe2)' EkBTR2
dd R
- -YD-
(20)
If hydrodynamic interaction is ignored, the last term on the right-hand side of Eq 20 is not present. For the simulations carried out in this work U(R > b ) is very small and also the second term on the right-hand side of eq 20 is negligible. Equation 16 then simplifies to
Diffusion-Controlled Target-Dimer Reactions
The Journal of Physical Chemistry, Vol. 88, No. 25, 1984
A
RS3.05
6155
I
which becomes equal to 47bDMin the absence of hydrodynamic interaction.
4. Details bf the Simulation In a given simulation, between 5000 and 10000 separate trajectories were calculated. Standard deviations were computed from the results for successive groups of 1000 trajectories. To begin a trajectory, the dimer was placed at R = b = 8 A and its initial orientation (ro) was selected at random from a Boltzmann distribution. The trajectory was terminated either when the reactive surfaces of the target and dimer collided or when R 1 q = 10 A. A reactive surface collision was considered to occur whenever R 5 3 A (3 8, equals the target sphere radius ( 2 A) plus twice the radius of the dimer subunits (0.5 A)) and 0 < 7r/2. When an unreactive surface collision occurred (0 > 7r/2), that particular dynamics step was ignored and a new step starting from the old coordinates was taken. In several simulations, both dimer subunits were considered reactive and the criteria for a reactive collision was simply that R S 3 A in those cases. As discussed previously:' the time step At must be sufficiently small so that displacements are (a) small relative to distances over which forces vary significantly, and (b) small relative to distances characterizing the domain of reflecting or absorbing boundary conditions. In the present work, CPU time was substantially reduced by varying At during a trajectory. In the simulations that included hydrodynamic interaction (HI), At was set equal to 2 fs for R < 5 %,or (1/160)D, = 6.5 fs for R > 5 A. The longer At corresponds to a net rms angular displacement of about 5" for the dimer. In simulations where HI was ignored, somewhat longer time ste s were taken for large R : AZ was set equal to 2 fs for R b. For b = 8 A, the longest time step corresponds to an rms angular displacement of about 36". In one simulation without HI, the smaller time steps employed in the H I simulations were used and the results were indistinguishable from the larger At run. Although a rmsdisplacement of 36" would produce large changes in the forces, the forces are weak at large R and contribute little to the net displacements. The reason the H I simulations are limited to shorter times is the restriction that d2 .- r R
