J . Phys. Chem. 1990, 94, 8794-8800
8794
of H 0 2 radicals in reactions 14 and 23. The overall homogeneous gas-phase production of HzS04takes on the order of 2-3 weeks. This means that the acid may be deposited some distance from the H2S source as the gas is transported in an air mass. Although some uncertainties remain, especially concerning the products of the HSO reactions, a general picture has been formed for the overall atmospheric H2S oxidation process.
Acknowledgment. This work was supported by NOAA as part of the National Acid Precipitation Assessment Program. We thank our colleague Dr. A. R. Ravishankara for helpful comments on this paper. Registry No. HS, 13940-21-1; HSO, 62470-71-7; O,, 10028-15-6; D,, 7782-39-0.
Kinetics of Diffusion-Influenced Reactions Studied by Brownian Dynamics Huan-Xiang Zhou Laboratory of Chemical Physics, National Institute of Diabetes and Digestive and Kidney Diseases, National Institutes of Health, Bethesda, Maryland 20892 (Received: February 22, 1990; In Final Form: May 30, 1990)
An algorithm is developed for calculating the time-dependent rate coefficient k ( t ) of diffusion-influenced reactions. This algorithm exploits the well-known relationship between the pair distribution function of the reactants, which defines k(r), and the survival probability, which can be easily simulated by using Brownian dynamics. Various boundary conditions are discussed. The efficiency of the algorithm is assessed. For several cases where k ( t ) is analytically known, this algorithm produces results that are in excellent agreement with their analytical counterparts over the entire time range. For other cases, the implementation of the algorithm does not depend on the particular geometry or interaction potential.
I. Introduction
has an equilibrium (i.e., Boltzmann) distribution around A
The kinetics .of diffusion-influenced reactions, characterized by the time-dependent rate coefficient k ( t ) , is analytically solvable for a limited number of cases.l-12 The rate coefficient k ( r ) is directly related to observed quantities in experiments such as ligand binding to macromolecule^,^^ fluorescence quenching,' ' - I 4 and the self-assembly of multisubunit proteins. In this paper, we present a computational algorithm for calculating k ( t ) based on Brownian dynamics simulations. The problem under consideration can be formulated as follows. We are interested in the diffusion-influenced reaction between two species, A and B. The concentration, C ( r ) ,of species B at r relative to species A normalized by C( m ) is defined as the pair distribution function, p ( r , r ) ,between the two species. The pair distribution function p(r,r) is assumed to satisfy the Smoluchowski equationi.iO,i la
p(r,O) = e-BU(I) (1.2) At t = 0, A and B are turned into reactive species (e.g., by photolysis). This is characterized by a radiation boundary condition2
d p ( r , t ) / d t = -V.J(r,t)
= V.De-BU(r)Ve~u(r)p(r,t)( 1 . I )
where 8 = (ken-', D is the relative diffusion coefficient of A and B (i.e., D = D, + D e ) , and U ( r ) is the potential of mean force between A and B. U p to time t = 0, A and B are inert and B ( I ) Smoluchowski, M. Z . Phys. Chem. 1917, 92, 129. (2) Collins, F. C.; Kimball, G.E. J . Colloid Sci. 1949, 4 , 425. (3) Weller. A. Z . Phys. Chem. 1957. 79, 551. (4) Carslaw, H. S.;Jaeger, J. C . Conduction o j H e a f in Solids, 2nd ed.; Oxford University Press: New York, 1959. ( 5 ) Hong. K. M.: Noolandi. J . J . Chem. Phys. 1978. 68, 5163, 5172. (6) Rice, S. A,; Butler, P.R.; Pilling, M. J.; Baird, J . K. J . Chem. Phys. 1979. 70, 4001. (7) (a) Pedersen, J. B. J . Chem. Phys. 1980, 72, 3904. (b) Pedersen, J . B.; Sibani, P. Ibid. 1981. 75, 5368. (c) Sibani, P.; Pedersen, J . B. Phys. Reo. Letr. 1983, 5 1 . 148. (8) Flannerq, M. R. Phys. Reo. 1982, ,425, 3403. (9) (a) Shoup. D.: Szabo. A. J . Elecrroanal. Chem. 1982, 140, 237. (b) Szabo, A.; Shoup, D.: Northrup, S. H.; McCammon. J . A . J . Chem. Phys. 1982, 7 7 , 4484. (IO) Rice, S. A. Diffusion-Limited Reactions. In Bamford, C. H.; Tipper, C. F. H.; Compton, R. G., Eds. Comprehcnsioe Chemical Kinetics; Elsevier: Amsterdam, 1985; Vol. 25. ( 1 1 ) (a) Szabo. A. J . Phys. Chem. 1989, 93,6929. (b) Agmon, N.: Szabo, A . J . Chem. Phys. 1990, 92, 5270. (12) Traytak, S. D. Chem. Phys. 1990, 140, 281. (13) The binding of carbon monoxide to protoheme provides an outstanding example. (14) Zhou. H. X . ; Szabo. A . J . Chem. Phys. 1990. 92, 3874.
-n.J(r,t)
= n.De-BU(r)VeBU(r)p(r,r) = K(r) p ( r , t ) ,
r E
r
(1.3a)
where r is the boundary between A and B, n is the normal of r, and K(r) is the space-dependence intrinsic bimolecular rate constant. Usually, part of r is reflecting; this is described by K(r) = 0. The absorbing extreme is described by K(r) = m . The outer boundary condition is given by p(m,t) = 1 (1.3b) The quantity of interest is the time-dependent rate coefficient defined by the equation
k ( t ) = - J - d S n.J(r,t) (1.4) r The reason that k ( t ) of eq 1.4 has the character of a rate coefficient is that the total concentrations of A, A(r), satisfies the kinetic equationi0>ll a dA(t)/dt = - k ( t ) C ( m ) A ( t ) (1.5)
in which it is assumed that B is in excess so that the change of C ( m ) with time is negligible. In principle, one can solve eq 1 . 1 for p(r,t) subject to the initial condition 1.2 and the boundary conditions 1.3 and integrate (eq 1.4) to obtain k ( t ) . This becomes formidable when the geometry or potential gets complicated. The purpose of this paper is to find an algorithm for calculating k ( t ) when arbitrary geometry and potential are present. In essence, our algorithm exploits the well-known relationship between the pair distribution function p ( r , t ) and the survival probability S(tlro)i5si6 Technically, S(tlr,) is the probability that reactants started (at time t = 0) at ro will survive (i.e., not react) at time 1. From the (15) Onsager, L. Phys. Reu. 1938, 54, 554. (16) Tachiya. M . J . Chem. Phys. 1978. 69, 2375
This article not subject to U S Copyright. Published 1990 by the American Chemical Societv
The Journal of Physical Chemistry, Vol. 94, No. 25, 1990 8195
Kinetics of Diffusion-Influenced Reactions simulation point of view, it is the fraction of reactants started at ro that are still alive at time t , which is easily calculated by using Brownian dynamics. For the steady-state rate coefficient, k ( m ) , Northrup, Allison, and McCammon" developed an elegant algorithm. This algorithm was recently rederived by using the relationship between the pair distribution function and the survival probability (eq 1.6).'* Besides the extension from the steady state to the transient regime, it will be seen that our algorithm also avoids some of the practical difficulties of the Northrup, Allison, and McCammon algorithm. The outline of this paper is as follows. In section 11, we develop the algorithm for calculating k ( t ) . In particular, we will discuss various boundary conditions, which present difficulties for Brownian dynamics simulations. To illustrate this algorithm, we apply it in section 111 to several cases. Finally, in section IV, we make some concluding remarks and point out future work. 11. The Algorithm We can write k ( t ) of eq 1.4 in terms of S(tlro), using eqs 1.3a and 1.6
(2.lb)
Figure 1 . T h e boundary
r.
This condition allows one to derive the relationship between the pair distribution function p(r,r) and the survival probability S(tlro), eq 1.6. Away from any boundary, if the diffusion coefficient D and F = -VU(r) are constant, the solution of eq 2.4 is
This suggests the Brownian dynamics algorithm2'
r = ro + DA@F(r0) + (2DAf)'/*R
where At is the time step and R is a vector of Gaussian random numbers with the following properties for its three components:
Here we have used the definitions
Equation 2.lc is the essential part of our algorithm. If one starts particles on the boundary r according to the distribution K(ro) exp[-@U(ro)] and counts the fraction of particles that survive reaction at time t , then k ( t ) is found. The dynamics of a particle is dictated by the conditional probability o(r,tlro,O), the probability that the particle is at r a t time 1, given that it was at ro at t = 0. The initial condition of u(r,tlro,O) is given by
a(r,O1ro,O) = 6(r-ro)
(2.2)
The relationship between S(tlro) and u(r,tlro,O) is then
S(tlro) = S d V u(r,tlro,O)
(2.3)
The conditional probability o(r,tlro,O) satisfies the same equation as p(r,t)1°
da(r,tlro,O)/dr = V.De-Bu(r)Ve~u(r)u(r,tlro,O) (2.4) and the same inner boundary condition
n ~ D e - B U ~ r ~ V ~ U ~ r ~ a (=r , K(r) f ~ r oa(r,rlro,O), ,O) rE
r
(2.5a)
(2.9a)
( R i R j )= 6,
(2.9b)
However, near a boundary, care needs to be taken to ensure that the boundary condition is satisfied. The ideal situation would be that one can find the exact solution of eq 2.4 subject to the boundary condition 2.5a and generate displacements from this distribution. We are interested in an algorithm, the implementation of which does not depend on the particular geometry of the boundary. The following strategy is taken. First, the radiation boundary condition, eq 2.5a, is converted to a reflecting boundary condition by putting a sink term in eq 2.4. Then the reflecting boundary condition is realized by generating r according to eq 2.8 and making a new attempt when r crosses the boundary. A comment is in order for each of the above two steps. Using a sink term to model reaction in the context of diffusion was discussed by Wilemski and Fixman22and by Northrup and H y n e ~ .The ~~ sink term approach and the radiation boundary approach are equivalent under appropriate conditions."a-14q22.23 As for implementing the reflecting boundary condition by making new attempts, it has been a common practice and, as will be seen below, works quite well. In what follows, we elaborate on these two steps. Now, instead of the conditional probability u(r,tlro,O) given by eqs 2.2, 2.4, and 2.5, we are interested in a modified conditional probability, p(r,tlro,O), that is governed by 1
(2.5b)
which says that a particle started (at time t = 0) at ro will have no chance of reaching infinity at any finite time. The conditional probability o(r,tlro,O) also satisfies the detailed balance condi
a(r,tlro,O)e-~u(ro)= u(ro,tlr,O)e-@u(r)
( R i )= 0
ap(r,tlro,O)/dt =
but the outer boundary condition
a(m,tlro,O) = 0
(2.8)
(2.6)
( 1 7) (a) Northrup, S.H.; Allison, S.A,; McCammon, J . A . J . Chem. Phys. 1984, 80, 1517. (b) Allison, S.A.; Northrup, S . H.; McCammon. J . A . J . Chem. Phys. 1985, 83, 2894. (18) Zhou, H. X. J . Chem. Phys. 1990, 92, 3092. (19) Gardiner. C . W . Handbook of Stochastic Methods, 2nd ed.; Springer-Verlag: Berlin, 1985. (20) Morse, P. M.; Feshbach, H. Methods of Theoretical Physics; McGraw-Hill: New York, 1953; Chapter 7 .
V.De-BU(r)VeBU(r)p(r,t~ro,O) - -p(r.tlro,O). 70.1
r E A (2.10a)
= V.De-Bu(r)VeBU(r)p(r,t~ro,O), r outside A
(2. lob) (2.1 OC)
n.De-BU(r)Veau(r)p(r,tlro,O)= 0, r E
r
(2.1 I )
where A is the region outside r where the normal of r is shorter than t (see Figure 1). In eq 2 . 1 0 ~it is understood that at each point, r , , along the normal of r at r, K(r,) is equal to K(r). (21) (a) Ermark, D. L. J . Chem. Phys. 1975,62,4189,4197. (b) Ermark, D. L.; McCammon, J . A . J . Chem. Phys. 1978, 69, 1352. (22) Wilemski, G.:Fixman, M. J . Chem. Phys. 1973, 58, 4009. (23) Northrup, S . H.; Hynes, J . T.J . Slat. Phys. 1978, 18, 91.
8796
The Journal of Physical Chemistry, Vol. 94, No. 25, I990
Zhou
--
-
6
Equations 2.2 and 2.5b still hold for p(r,t)ro,O). When t 0, p(r,tlro,O) o(r,tlro,O). To see this, first note that when e 0, eq 2.10b reduces to eq 2.4. Then integrate eq 2.10a along the normal of r O = ,.De-su(r)Desu(r)p(r,rlro,O) -
( I / € ) x f d n l 4 r 1 ) p(rl,tlro,O),
r E
r'
(2.12)
i
where I" is the outer boundary of A, which is c away from r. Equation 2.1 1 has been used in deriving eq 2.12. Clearly, when c 0, eq 2.12 reduces to eq 2.5a. As a consequence, p(r,tlro,O) u(r,tlro,O). In practice, t 2 / D r a. When eqs 2.1 and 2.14 are then used of r (I”). This will be the situation in all the cases that we study. to calculate the steady-state rate coefficient, it is found that (eq We emphasize that this is just a matter of convenience; for the ,441 sake of starting particles one can put K(ro) into a form of potential. In any event, if in addition no force is present, then particles are Drc/Ka2)erc/a - I]-’ (2.29) & ( a ) = 47rDr,[( I Otherwise they have to be started started uniformly over rl (rl’). and (eqs AIO) according to exp[-PU(ro)]. In general this can be done through a Monte Carlo sampling, but specific cases may have simpler &,(a) = + (a + solutions (see below). We mention again that particles are terc ) - ~ ( D T ) - ’ / cosh ~u [ ( D ~ ) - l / ~+ c ] ( a + e)-’ sinh [ ( D T ) - ~ / ~ C ] ] minated after their lifetimes exceed some preset cutoff. (2.30a) We now comment on the start and stop mechanisms just described. They are directly related to the efficiency of the al(2.30b) = k ( m ) [ I - Y t 0(€2)] gorithm. In the Northrup, Allison, and McCammonI7 algorithm for calculating k ( m ) , particles are started on a spherical surface r = b, beyond which the potential U(r) is spherically symmetric. -, Then they are allowed to evolve until they react or touch a spherical surface r = q, beyond which the steady-state pair disih [ (D T ) - ’ / ~ ~ ] tribution function p(r,m) is spherically symmetric (see ref 18). As discussed in ref 18, this algorithm becomes cumbersome when (2.30~) a strongly noncentrosymmetric force is present or when the reactive portion is surrounded by a large inert part. Actually, a modified ( ~ / 1 2 D I/a r , / 4 ~ ~ ) e ? c /(~~-/ 1 2 D I/a + rc/a2) version of the Northrup, Allison, and McCammon algorithm was 1 used to deal with such situations.2s In our algorithm, particles (1 Drc/Ka2)~c/a- 1 (2.30d) S,(tlm) = 1
(2.28b)
-
-
+
+
L\-
+
+
‘ I
+
+
