J. Phys. Chem. B 2001, 105, 6017-6024
6017
An Efficient Brownian Dynamics Method for Evaluating Inertial Dynamic Effects on Diffusion-Influenced Reactions Seongeun Yang, Heekyung Han, and Sangyoub Lee* School of Chemistry and Molecular Engineering, and Center for Molecular Catalysis, Seoul National UniVersity, Seoul 151-747, South Korea ReceiVed: January 22, 2001; In Final Form: March 27, 2001
We present a very efficient Brownian dynamics (BD) method for evaluating inertial dynamic effects on diffusion-influenced reaction kinetics. The method requires only about 10 000 trajectories to yield a reliable estimate of the time-dependent rate coefficient. Another key advantage of the present BD method is that it requires only one set of BD simulations to calculate the rate coefficient for any value of the intrinsic rate constant. By using the proposed method, we show that the subtle inertial dynamic effect that may arise from the coupling between the center-of-mass and the relative motions of reactants is negligible.
1. Introduction The theoretical description of the chemical reaction rate influenced by the diffusive encounter rate of reactant molecules has usually been based on the reaction-diffusion equation.1-6 Various factors that may affect the diffusion-influenced reaction kinetics have been studied, which include intermolecular forces among the reactants and solvent caging, hydrodynamic repulsion, anisotropic reactivity and rotational diffusion, gating mode, higher-order concentration effects, and so on.2-10 Although the diffusion equation (DE) description has provided a better explanation of experimental observations and computer simulation results by including these factors, it still fails to give an exact account of the short-time scale reaction dynamics. In an attempt to assess the validity of the DE description, Dong, Baros, and Andre11 (DBA) carried out a molecular dynamics (MD) simulation study on the kinetics of fluorescence quenching (A + B f B, with A the excited fluorophor and B the quenchers). They considered a model reaction system that consists of a single molecule of A and many molecules of B in a bath of solvent molecules. All the molecules were modeled as hard spheres with identical mass and size, and some of them were tagged as reactant molecules. Starting from 500 identical hard spheres at equilibrium, they evaluated the time-dependent survival probability of A for differing numbers of B molecules. The results were then compared with the prediction of the simplest version of the Smoluchowski theory, in which the reactants react infinitely fast at contact and the potential of mean force and the hydrodynamic interaction between them are neglected. Major discrepancies were found. The simulated survival probability decays much faster than the Smoluchowski theory predicts. DBA attributed this discrepancy to the nondiffusional dynamic behavior of reactant molecules over short time scales.11,12 However, the simplest Smoluchowski theory involves many other assumptions that can be invalidated in interpreting the DBA simulations. Indeed, Zhou and Szabo showed that the discrepancy between the theory and simulations can be surprisingly small except for a high concentration of B molecules if the theory employs the radiation boundary condition at contact with an intrinsic rate constant determined by the low-density
collision frequency and the potential of mean force determined by the pair distribution function.13 It is evident that a separate assessment of the validity of every assumption made in the Smoluchowski theory is desirable to pinpoint the major cause of the disagreement between the theory and MD simulations. In an attempt along these lines, we recently investigated the effect of excluded volumes among like reactant molecules.14 By comparing the results of an analytical theory based on the DE with those obtained from diffusive-regime Brownian dynamics (BD) simulations, we have found that the excluded-volume effect is the major cause of the discrepancy between the improved version of the Smoluchowski theory and MD simulations as observed by Zhou and Szabo for a high concentration of B molecules. The purpose of the present work is to investigate the effect of nondiffusional motions of reactant molecules on the reaction kinetics. The nondiffusional dynamic behavior should be analyzed from two aspects. The first is the non-Markovian dynamic effect due to finite relaxation times of the solvent motions, and the second is the inertial dynamic effect due to the finite velocity relaxation time of the reactant molecules. These aspects have been investigated by using the generalized diffusion equation12,15 and the ordinary or generalized FokkerPlanck equation (FPE).15-19 The non-Markovian effect has been found to be negligible at the usual fluid density.15 In this study, we therefore focus on the inertial dynamic effect within the theoretical framework based on the ordinary FPE. Unfortunately, analytical solutions to the FPE are available only for a few simple cases, and nontrivial problems involving reactive (absorbing or radiative) boundary conditions cannot be solved exactly with presently known mathematical techniques.20,21 Harris16 derived an approximate analytic expression for the Laplace transform of the time-dependent rate coefficient for the diffusion-controlled reaction in three dimensions by considering the capture rate of Brownian particles by a completely absorbing sphere located at the origin of coordinates. He employed the lowest-order approximation in a general halfrange moment expansion method, which takes the phase-space distribution function as a product of the bimodal radial density distribution and the Maxwell velocity distribution at all times. The dynamic correlation between the position and velocity
10.1021/jp0102419 CCC: $20.00 © 2001 American Chemical Society Published on Web 06/05/2001
6018 J. Phys. Chem. B, Vol. 105, No. 25, 2001 spaces is taken into account through the bimodality of the radial density distribution. However, the assumption that the velocity distribution is Maxwellian may not hold in the kinetic boundary layer near the absorbing surface. Since the fate of the reactants is determined in this region within a very short time, the resulting error can be significant. Therefore, a more accurate evaluation of the inertial dynamic effects based on appropriate computer simulation methods is desirable. When the concentrations of reactant molecules are low, their reaction dynamics can be described by those of an ensemble of independent reactant pairs. Then, one usually considers only the relative motion of a reactant pair to determine the reaction kinetics. This approach is valid as long as the center-of-mass (CM) motion and the relative motion are decoupled. For example, if the thermal motions of reactants can be described by the Smoluchowski equation, the CM motion does not affect the relative motion. However, when the dynamics needs to be described by the Fokker-Planck equation and the reactants do not have identical masses and friction coefficients, the CM motion couples with the relative motion.22 This coupling is another kind of subtle inertial dynamic effect, but its influence on the reaction kinetics has not been investigated before. In the present work, we will evaluate the inertial dynamic effect by carrying out BD simulations for a system of two spherical reactants. We have recently developed a very efficient BD method for calculating the time-dependent rate coefficient, which is applicable when the dynamics of the reactants is governed by the DE.23 In section 2, we extend the method to deal with the case when the reactant dynamics is governed by the FPE. The method can be easily generalized to reactions of more complicated types involving reactants of arbitrary shape and interaction potential. Details of numerical implementation are described in section 3. In section 4, we evaluate the inertial dynamic effect on the reaction between identical molecules of the same mass, size, and diffusion constant. We calculate the variation of the rate coefficient when the mass of the reactants changes but their size and diffusion constant are fixed. Since the DE description gives the same rate coefficient in this case independent of the reactant mass, the result will provide a direct measure of the inertial dynamic effect. In section 5, we evaluate the suspected effect of coupling between the CM and the relative motions on the reaction between unlike reactants. It is found that this inertial coupling effect on the reaction rate is practically negligible. 2. The BD Method We will consider the irreversible bimolecular reaction, A + B f product(s), in the low reactant concentration limit. We first generalize the Lee-Kaplus (LK) formalism24 for treating the diffusion-influenced reaction kinetics to include the inertial dynamic effect, and derive an expression for the timedependent rate coefficient kf(t) in terms of the returning probability (RP), which is the probability that a pair of reactant molecules in a reaction zone will be found again in the zone at a later time under the condition that the reaction is absent. The RP is the key dynamic quantity that is evaluated from BD simulations. We denote the phase-space coordinates (position and velocity) of molecules A and B by X ) (rA,vA) and Y) (rB,vB), respectively. Let ΦA(X,t) be the one-particle reduced distribution function (RDF) that gives the number density of A molecules at rA with the velocity given by vA. ΦB(Y,t) is defined similarly, and ΦAB(X,Y,t) is the two-particle RDF representing the product of number densities of A and B molecules at rA and rB with
Yang et al. velocities vA and vB. The evolution equations governing these RDF’s are24
∂ Φ (X,t) ) LAΦA(X,t) ∂t A
∫ dY SAB(X,Y)ΦAB(X,Y,t)
∂ Φ (Y,t) ) LBΦB(Y,t) ∂t B
∫ dX SAB(X,Y)ΦAB(X,Y,t)
(2.1)
(2.2)
∂ Φ (X,Y,t) ) LABΦAB(X,Y,t) - SAB(X,Y)ΦAB(X,Y,t) ∂t AB
∫ dY′ SAB(X,Y′)ΦABB(X,Y,Y′,t) -∫ dX′ SAB(X′,Y)ΦABA(X,Y,X′,t)
-
(2.3)
LA, LB, and LAB are the Fokker-Planck operators governing the nonreactive dynamics of A and B molecules,
LR ) -vR ‚
(
∂ ∂ + γR ‚v + ∂rR ∂vR R γRkBT ∂ ∂ ‚ mR ∂vR ∂vR
LAB ) -vA ‚
(R ) A,B) (2.4)
∂ 1 ∂UAB ∂ ∂ + γA ‚ vA + ‚ + ∂rA ∂vA mA ∂rA ∂vA
) (
γAkBT ∂ ∂ ∂ ∂ ‚ + -vB ‚ + γB ‚v + mA ∂vA ∂vA ∂rB ∂vB B γBkBT ∂ 1 ∂UAB ∂ ∂ ‚ + ‚ (2.5) mB ∂rB ∂vB mB ∂vB ∂vB
)
where γR and mR are the friction coefficient and the mass of molecule R, kBT is the product of the Boltzmann constant and the absolute temperature, and UAB(X,Y) is the potential of mean force between A and B molecules. SAB(X,Y) is the sink function that describes the depopulation rate of the reactant molecules due to reaction at the configuration (X,Y). The last two terms in eq 2.3, involving the three-particle RDFs ΦABB and ΦABA, arise from the competitive reactions of surrounding reactants with a given pair of A and B molecules. To obtain a manageable set of kinetic equations, we will employ superposition approximation as given in eqs 2.9 and 2.10. We assume that reactant molecules are distributed in equilibrium at t ) 0. Then, the one-particle RDFs remain uniform over the spatial coordinates at any t (g0), since in the absence of external fields reaction occurs homogeneously throughout the medium and the microscopic environment of any reactant molecule is isotropic on the average. We can thus write
ΦA(X,t) ) [A]PvA(vA,t)
(2.6)
ΦB(Y,t) ) [B]PvB(vB,t)
(2.7)
where [A] and [B] are the bulk number densities of A and B at time t. PvR(vR,t) is the probability density that a molecule of species R which has not reacted until t has the velocity vR and is assumed to be given by an equilibrium distribution, i.e., the Maxwellian distribution at t ) 0. If the reactivity depends on vR, PvR(vR,t) will deviate from the equilibrium value for t > 0 since there is a preferential removal of R molecules with velocities in some particular range due to the reaction.
An Efficient Brownian Dynamics Method
J. Phys. Chem. B, Vol. 105, No. 25, 2001 6019
The two-particle RDF ΦAB is related to the pair correlation function FAB by
ΦAB(X,Y,t) ) [A][B]PvA(vA,t)PvB(vB,t)FAB(rBA,vA,vB,t) (2.8) where rBA ) rB - rA. We will denote the set of variables (rBA,vA,vB,t) collectively by Λ. Employing the reactive superposition approximation,24,25 the three-particle RDFs ΦABB and ΦABA are expressed as
ΦABB(X,Y,Y′,t) ) [A][B]2PvA(vA,t)PvB(vB,t)PvB(v′B,t)FAB(Λ,t)FAB(Λ′,t) (2.9) ΦABA(X,Y,X′,t) )
and from eq 2.12, we have
∂ ([A][B]fAB) ) -kf(t)[A][B]2fAB - kf(t)[A]2[B]fAB + ∂t ∂fAB [A][B] (2.16) ∂t Combining eqs 2.15 and 2.16, we obtain
∂fAB ) LABfAB(Λ,t) - SAB(Λ)fAB(Λ,t) + ∂t
∫ dr′B dv′B SAB(Λ′)PvB (v′B,t)FAB(Λ′,t)]fAB(Λ,t) + [A][kf(t) - ∫ dr′A dv′A SAB(Λ′′)PvA (v′A,t)FAB(Λ′′,t)]fAB(Λ,t)
[B][kf(t) -
[A]2[B]PvA(vA,t)PvB(vB,t)PvA(v′A,t)FAB(Λ,t)FAB(Λ′′,t) (2.10)
(2.17)
where Λ′ ≡ (r′B - rA,vA,v′B) and Λ′′ ≡ (rB - r′A,v′A,vB). The reactive superposition approximation neglects the correlation between like particles. As shown by Kotomin and Kuzovkov,6 like particles tend to cluster as the bimolecular reaction proceeds, since a particle located in an environment where nonreacting like particles are more densely populated than reacting partners would have a better chance of surviving. This kind of likeparticle correlation effect can be important in the long-time limit. However, the more important aspect of like-particle correlations is the effect of excluded volumes;14 the strong repulsive interactions between like reactant particles prevent them to overlap. This excluded-volume effect will be important when the reactant concentration is high. Also, it is known that the superposition approximation is less accurate in lower dimensional space.6 Substituting eqs 2.6 and 2.8 into eq 2.1 and integrating the resulting equation over vA, we obtain
The last two terms on the right-hand side of eq 2.17 almost vanish when the inertial coupling between the CM and the relative motions can be neglected, since in such case FAB(Λ,t) as well as the sink function is a function of the relative position rBA and the relative velocity vBA() vB - vA) only. In any case, these are the higher-order concentration terms and will be neglected hereafter. That is, we have in the low reactant concentration limit
∂ {[A] ∂t
∫ dvAPvA} ) [A]∫ dvALAPvA [A][B]∫ dΛ SAB(Λ)PvA PvBFAB(Λ,t)
(2.11)
where we have noted that SAB(X,Y) ) SAB(Λ). Since ∫ dvAPvA ) 1, and ∫ dvALAPvA ) 0, we obtain the rate law from eq 2.11:
d[A] ) -kf(t)[A][B] dt
∫ dΛ SAB(Λ)fAB(Λ,t)
(2.18)
To go further, we now consider a sink function of the form
SAB(Λ) )
{
kR if rBA‚vBA < 0 and σ e |rBA| e R (2.19) 0 otherwise
kR is an intrinsic rate parameter, σ is the contact distance between A and B molecules, and R is the radius of the outer boundary of the reaction zone. In the limit of kR f ∞ and R f σ, the problem of solving eq 2.18 is equivalent to that of solving the equation ∂fAB/∂t ) LABfAB(Λ,t) with the absorption boundary condition considered by Harris16 and Naqvi et al.17 By applying the operator algebra, we can obtain a formal expression for the time-dependent rate coefficient as24
kf(t) ) k°f - (k°f)2
(2.12)
∫0t dτ∆(τ)
(2.20)
Here, k°f is the rate constant that would be observed if the pair correlation between A and B molecules is maintained at equilibrium,
The rate coefficient kf(t) is given by
kf(t) )
∂fAB ) LABfAB(Λ,t) - SAB(Λ)fAB(Λ,t) ∂t
(2.13)
with
k°f )
∫ dΛ SAB(Λ)gAB(Λ)
(2.21)
with
fAB(Λ,t) ≡
PvA(vA,t)PvB(vB,t)FAB(Λ,t)
(2.14)
To evaluate kf(t) from eq 2.13, we need to determine fAB(Λ,t). An evolution equation for fAB can be obtained from eq 2.3. Substituting eqs 2.8-2.10 into eq 2.3, we obtain
∂ ([A][B]fAB) ) [A][B]LABfAB(Λ,t) ∂t [A][B]SAB(Λ)fAB(Λ,t) -
∫ dr′B dv′B SAB(Λ′)PvB (v′B,t)FAB(Λ′,t)fAB(Λ,t) [A]2[B]∫ dr′A dv′A SAB(Λ′′)PvA (v′A,t)FAB(Λ′′,t)fAB(Λ,t)
[A][B]2
(2.15)
gAB(Λ) ) e-UAB(r)/kBT φA(vA)φB(vB) φR(vR) )
( ) ( ) mR 2π kBT
mRvR2 (R ) A, B) 2kBT
3/2
exp -
(2.22)
The memory kernel ∆(τ) is the key dynamic function. For the sink function given by eq 2.19, the following approximate relation can be obtained:
ˆ 0(z)] ∆ ˆ (z) = ∆ ˆ 0(z)/[1 + k°f ∆
(2.23)
We denote the Laplace transform of a function f(t) by ˆf(z). In eq 2.23, ∆0(t) is the memory kernel in the absence of reaction
6020 J. Phys. Chem. B, Vol. 105, No. 25, 2001
Yang et al.
that is given by
The position of the particle is then given by
Vrx∆0(t) )
∫rx dΛ 〈G(Λ,t|Λ0)〉0
(2.24)
Here, ∫rx dΛ denotes an integration over the reaction zone that is represented by the condition rBA‚vBA < 0 and σ e |rBA| e R. Vrx is a quantity that may be roughly described as the “volume” of the reaction zone in the Λ space:
Vrx )
∫rx dΛ gΑΒ(Λ) = 2πσ2∆r exp[-UAB(σ + 21∆r)/kBT]
r(t + ∆t) ) r(t) +
G(Λ,t|Λ0) ) etLAB(Λ)δ(Λ - Λ0)
(2.26)
where LAB(Λ) is the evolution operator governing the thermal motion of the reactant pair in the Λ space (see eq 5.1). In eq 2.24, the angular bracket denotes a Boltzmann weighted average over the initial distribution within a reaction zone:
∫rx dΛ0 gAB(Λ0)(‚‚‚) 〈(‚‚‚)〉0 ) ∫rx dΛ0 gAB(Λ0)
〈B22〉 )
3. Numerical Implementation The overall procedure for generating and analyzing the BD trajectories follows the same lines as described in the previous work.23 Trajectories for a pair of reactant molecules A and B are initiated in a reaction zone satisfying the condition rBA‚vBA < 0 and σ e |rBA| e R. First, with the molecule A located at the origin of coordinates, the initial position of the molecule B is sampled on the x-axis according to the distribution r2e-UAB(r)/kBT/ ∫Rσ dr r2 e-UAB(r)/kBT (r ) |rBA|). Then, the velocities of the molecules are sampled according to the Maxwellian distributions φA(vA) and φB(vB). If the velocities satisfy the condition rBA‚ vBA < 0, they are accepted. Otherwise, the sampling is repeated. In this way, the averaging operation in eq 2.27 is automatically conducted. The motion of the reactant particles obeys the Langevin equation, and we use the BD move algorithm developed by Ermak and Buckholz.26 For a particle with mass m and friction coefficient γ, the velocity after a time step ∆t is given by
v(t + ∆t) ) v(t)e-γ∆t +
F(t) (1 - e-γ∆t) + B1 mγ
(3.1)
F(t) is the systematic force exerted on the particle, and B1 is the random velocity change that is chosen from a Gaussian distribution with a mean value of zero and a variance of
〈B21〉 )
3kBT (1 - e-2γ∆t) m
(3.2)
[
]
1 - e-γ∆t γ∆t 2 mγ2 1 + e-γ∆t
6kBT
(3.4)
The random vectors B1 and B2 are generated by the IMSL subroutine DRNNOA.27 The time step size ∆t must be small so that the change in the systematic force during ∆t can be neglected. Also, for an improved evaluation of the inertial effect, ∆t must be smaller than the velocity relaxation time. In addition, to obtain an accurate estimate of the returning probability at short times, the trajectories started from the middle part of the reaction zone should not leave the reaction zone in just one or two steps. The initial time step size (∆t)0 is determined by the relation,
(∆t)0 )
(2.27)
Therefore, the quantity on the right-hand side of eq 2.24 represents the returning probability that a pair of reactant molecules located in the reaction zone at t ) 0 will be found again in the reaction zone at a later time t under the condition that the reaction is absent. We can calculate this returning probability from BD trajectories and thereby the time-dependent rate coefficient kf(t).
)
where the random displacement B2 is chosen from a Gaussian distribution with a mean value of zero and a variance of
(2.25) We are assuming that the width of the reaction zone, ∆r(≡R σ), is very small. It can be noted that k°f ) kRVrx. G(Λ,t|Λ0) is the Green’s function in the absence of reaction:
(
1 v(t + ∆t) + v(t) γ 2F(t) 1 - e-γ∆t F(t) ‚ + ∆t + B2 (3.3) mγ 1 + e-γ∆t mγ
1 (∆r) C 2D
2
(3.5)
where D is the relative diffusion constant. The value of C should preferably be larger than 10. Away from the reaction zone (r > rb), the time step size is varied according to
(
1 (r - rb) ∆t ) (∆t)0 + Int C ′ 2D
)
2
(3.6)
where Int(x) gives the largest integer value smaller than x. The constant C ′ and the separation rb should be determined such that the trajectory cannot return to the reaction zone in a single step from a distance larger than rb. We assume that the reactant molecules are at the instant of collision, when |σ - |rBA|| is smaller than a preset value, 10-4σ. Then, the velocities of reactant molecules are modified according to28
v′A ) vA +
2mB (rˆ c ‚v )rˆ c and mA + mB BA BA BA 2mA (rˆ c ‚v )rˆ c (3.7) v′B ) vB mA + mB BA BA BA
where rˆ cBA ) rcBA/|rcBA| and rcBA is the distance vector at the collision. On the other hand, if the reactant molecules overlap too strongly in a given step so that the value of σ - |rBA| is larger than the preset value, the step is withdrawn and retried in two steps. First, the reactant molecules retry their movements from the old positions and velocities with a smaller time step size ∆t1, so that the overlap between the reactants becomes smaller than the preset value. Then, the next movement is searched with the time step size given by ∆t - ∆t1. Trajectories are terminated when the time length exceeds some cutoff. The cutoff time T should be sufficiently large that the returning probability diminishes to almost zero for t > T. In practice, this cutoff time T can be determined from a preliminary calculation. From the record of N trajectories of time length T, the returning probability can be calculated as
An Efficient Brownian Dynamics Method
J. Phys. Chem. B, Vol. 105, No. 25, 2001 6021
Pret(t) ≡ Vrx∆0(t) )
1 (number of trajectories found in the reaction N zone at time t) (3.8)
The reaction-free memory kernel ∆0(t) so obtained is Laplacetransformed numerically by using the IMSL subroutine DQ2AG.27 However, ∆0(t) calculated from a finite number of trajectories may not be smooth because of the statistical errors, especially at long times when the returning probability is far less than unity. Hence it needs to be smoothed out first. This could be done by averaging the returning probability, obtained from eq 3.8, over a time interval. In addition, to obtain a value of the returning probability function at an arbitrary time, interpolation of the discrete data points is necessary. We use the IMSL subroutines DCSSED and DCSINT for data smoothing and interpolation.27 Finally, the time-dependent rate coefficient kf(t) is obtained by the inverse Laplace transformation from
ˆ 0(z)] kˆ f(z) ) z-1k°f/[1 + k°f∆
(3.9)
We use the IMSL subroutine DINLAP for the numerical inverse Laplace transformation.27
Figure 1. Comparison of the BD results (dashed curves) with the exact ones (solid curves) calculated from eq 4.2, for the four cases with ∆r ) 0.1, 1.0, 2.0, and 5.0 Å. Since the exact result is available only for the high friction case, the value of the friction coefficient γR is set to 108 ns-1. The upper curves are for the case with f ) k °f /[4πσ(DA + DB)] ) 1 and the lower curves are for the case with f ) 5.
4. Reactions between Identical Hard Spheres In this section, we examine the effect of inertial dynamics on the reaction between identical hard spheres. We take the reaction occurring in the high friction limit (with large γ) as the reference for measuring the magnitude of the inertial effect. Szabo29 considered a reaction model in which the motions of reactants are governed by the DE and the reaction sink function has the form,
SAB(|rBA|) )
{
kR if σ e |rBA| e R 0 otherwise
(4.1)
In the high friction limit our reaction model reduces to this one, except that the equilibrium rate constant k°f for our model is one-half of that for Szabo’s model, since in our model the reaction can occur only when the reactants head for each other. The rate coefficient for Szabo’s model is given by29
kˆ f(z) ) Vrx
( ) ( )
kR kR 2 K(z)kˆ S(z,R) + (4.2) kR + z kR + z K(z) + zkˆ S(z,R)
where
kˆ S(z,R) ) (4πDR/z)[1 + (zR2/D)1/2] K(z) ) 4πDR
[
]
λ∆r cosh λ - (∆r - λ2σR/∆r)sinh λ λσ cosh λ + ∆r sinh λ
with λ ) ∆r[(kR + z)/D]1/2. The time-dependent rate coefficient can be obtained from eq 4.2 by the inverse Laplace transformation. The present BD method is based on an approximate relation between the memory kernel in the presence of reaction and that in the absence of reaction given by eq 2.23. This relation is exact for the δ-function-like sink function, but is expected to be accurate only for the case with a narrow reaction zone. Therefore, we first examine the accuracy of the present BD method with respect to the reaction zone width, ∆r. We consider the reaction between two identical hard spheres of 5 Å radius
with a diffusion constant of 25 Å2/ns, and neglect the potential of mean force. Since an exact result is available only for the high friction case, we set the value of the friction coefficient γR to 108 ns-1. In Figure 1, BD results (dashed curves) for the four cases with ∆r ) 0.1, 1.0, 2.0, and 5.0 Å are compared with the exact ones (solid curves) calculated from eq 4.2. One can see that the present BD method is sufficiently accurate even when ∆r ) 5.0 Å. In each of the four figures, the upper curves are for the cases with k°f ()kRVrx) having the same magnitude as 4πσD (D ) DA + DB), while the lower curves are for the cases with k°f five times larger than 4πσD. In the BD simulations, the initial time step size (∆t)0 was set by using eq 3.5 with C ) 100. Then, the time step size was varied according to eq 3.6; the value of C ′ was set to 100, and rb was set to (R + 1) Å (for ∆r ) 0.1 and 1.0 Å) or (R + 2) Å (for ∆r ) 2.0 and 5.0 Å). In each case, we have been able to obtain statistically reliable results from about 10 000 trajectories. We then estimate the effect of inertial dynamics by using the same reaction model: the reaction between two identical hard spheres of 5 Å radius with a diffusion constant of 25 Å2/ ns. The potential of mean force is neglected again, and the width of the reaction zone, ∆r, is taken to be 0.1 Å. The initial time step size (∆t)0 is set to 1.0 fs, as calculated from eq 3.5 with C ) 100. Then for r > rb ) 11 Å, the time step size is varied according to eq 3.6 with C′ ) 100. To examine the inertial effect, we have fixed the values of the reactant radius and the diffusion constant but have varied the mass and the friction coefficient according to the Stokes-Einstein relation, mRγR ) kBT/DR (a ) A, B), with T ) 300 K. In all cases, we have been able to obtain reliable results from about 10 000 trajectories. Figure 2 displays the inertial dynamic effect on the returning probability. The solid curve represents the exact result known for the diffusion limit (γR f ∞).23 For the reaction sink function given by eq 2.19, the initial velocity distribution satisfies the condition, rBA‚vBA < 0. However, in the diffusion limit, it immediately relaxes to the Maxwellian distribution. Hence the returning probability becomes one-half at t ) 0+, and its value for t > 0 can be calculated from eq 4.3 of ref 23 by simply
6022 J. Phys. Chem. B, Vol. 105, No. 25, 2001
Yang et al. coordinates. We have
LAB ) LCM + LRM + LC with
ΓkBT ∂ ∂ ∂ ∂ ‚ LCM ) -V‚ + Γ ‚V + ∂R ∂V M ∂V ∂V LRM ) -v ‚
Figure 2. Inertial dynamic effect on the returning probability Pret(t) for the reaction between identical hard sphere reactants. The diffusionlimit result is represented by the solid curve, while the BD simulation results are represented by the dot-dashed curve (γR ) 108 ns-1) and the dashed curve (γR ) 2 × 103 ns-1).
Figure 3. Inertial dynamic effect on the rate coefficient for the reaction between identical hard sphere reactants. The upper curves are for the case with f ) k°f/[4πσ(DA + DB)] ) 1 and the lower curves are for the case with f ) 5. The diffusion-limit results calculated from the analytic expression in eq 4.2 are represented by the solid curves, while the BD simulation results are represented by the dot-dashed curves (γR ) 108 ns-1) and the dashed curves (γR ) 2 × 103 ns-1).
dividing it by two. The dot-dashed curve is the BD simulation result for a high friction case with γR ) 108 ns-1. As expected, it almost coincides with the solid curve. However, the γR value of 108 ns-1 is rather unphysical. A reasonable value of the molecular mass that is appropriate for the size of the molecule under consideration is about 5000 g/mol. With mR ) 5000 g/mol, we have γR ) 2 × 103 ns-1. The dashed curve in Figure 2 represents this case. Its deviation from the solid curve thus provides a physically meaningful estimate of the inertial dynamic effect. We note that the inertial dynamics of the reactants influences the returning probability on the subpicosecond time scale. Figure 3 shows that the inertial dynamic effect on the rate coefficient persists for a long time. The upper three curves are for the cases with f ≡ k°f/4πσD ) 1, while the lower three curves are for the cases with f ) 5. The solid curves represent the time-dependent rate coefficients calculated from eq 4.2. The dot-dashed curves are the BD simulation results for the high friction case with γR ) 108 ns-1, while the dashed curves represent the results for the case with γR ) 2 × 103 ns-1. The reaction rate decreases as the mass of the reactant increases while keeping the product of its size and diffusion constant fixed. This trend is in accord with the prediction of Harris.16 5. Reactions between Unlike Reactants The Fokker-Planck operator LAB, which was given by eq 2.5, can be rewritten in terms of the CM and the relative
γkBT ∂2 ∂ 1 ∂UAB ∂ ∂ +γ ‚v+ ‚ + ∂r ∂v µ ∂r ∂v µ ∂V2
LC ) (γB - γA)
[
]
mAmB 2kBT ∂ ∂ ∂ ∂ v‚ +V‚ + ‚ 2 ∂V ∂v M ∂V ∂v M (5.1)
where (R,V) are the CM coordinate and velocity and (r,v) are the relative coordinate and velocity (r ≡ rB - rA). M ) mA + mB and Γ is the effective friction coefficient for the CM motion given by Γ ) (mAγA + mBγB)/M. µ is the reduced mass and γ is the reduced friction coefficient given by γ ) (mAγB + mBγA)/ M. LCM is the Fokker-Planck operator for the CM motion, while LRM is that for the relative motion. Owing to the presence of the term LC, these motions are in general coupled. When the two reactants have the same friction coefficient, the coupling vanishes. Rodger and Sceats have shown that the coupling vanishes also when the CM velocity distribution is symmetric.22 It has usually been assumed that the reaction dynamics between two hard spherical molecules can be described by considering the relative motion only.12,15-17 However, eq 5.1 shows that this may not be true for the case with unlike reactants. We now evaluate the suspected effect of coupling between the CM and the relative motions from BD simulations. As the reference reaction, we take the fixed target model, in which one reactant is fixed at the origin of coordinates and the other reactant with a reduced mass µ and friction coefficient,
γrel )
kBT µ(DA + DB)
(5.2)
moves around it. Note that as for the friction coefficient associated with the relative motion, γrel given by eq 5.2 is a better choice than γ as introduced in eq 5.1. The relative Fokker-Planck operator LRM is meaningless unless the coupling with the CM motion is properly considered. For example, in the high-friction regime the Fokker-Planck equation, ∂f(r,v,t)/ ∂t ) LRMf(r,v,t), reduces to a diffusion equation with the relative diffusion coefficient given by kBT/µγ. Unless the two particles are identical, this deviates from the correct one given by D ) DA + DB. Effects of the coupling of the CM motion with the relative motion are measured by comparing the results of the moving target model, in which both reactants move, with those obtained for the fixed target model. Figure 4 displays the results. For the four subplots, we have fixed the values of D and σ, but have varied the reactant mass ratio fm ) mB/mA with D ) 50 Å2/ns, σ ) 10 Å, and mA ) 1000 g/mol. fm ) 1.0 in panel a, 0.2 in panel b, 0.1 in panel c, and 0.01 in panel d. The diffusion constant DR and the radius σR of each individual molecule are calculated from the following relations:
rA ) σ/(1 + fm1/3), rB ) σ - rA, DA ) D/(1 + fm-1/3), and DB ) D - DA (5.3)
An Efficient Brownian Dynamics Method
J. Phys. Chem. B, Vol. 105, No. 25, 2001 6023 This indicates that one need not bother with the suspected effect of coupling between the CM and the relative motions in considering the reaction between unlike reactants. 6. Conclusions
Figure 4. Effect of the coupling between the center-of-mass and relative motions on the returning probability Pret(t). The reactant mass ratio has been varied; fm ) mB/mA ) 1 in panel a, 0.2 in panel b, 0.1 in panel c, and 0.01 in panel d. Values of other parameters used are described in the text. The solid curve represents the fixed target result, while the dashed curve represents the moving target result.
Figure 5. Effect of the coupling between the center-of-mass and relative motions on the reaction rate coefficient. The reactant mass ratio has been varied; fm ) mB/mA ) 1 in panel a, 0.2 in panel b, 0.1 in panel c, and 0.01 in panel d. The upper curves are for the case with f ) k°f/[4πσ(DA + DB)] ) 1 and the lower curves are for the case with f ) 5. Values of other parameters used are described in the text. The BD simulation results for the fixed target model (solid curves) can hardly be distinguished from those for the moving target model (dashed curves).
The friction coefficients are then calculated from the Einstein’s relation, γR ) kBT/mRDR. In the figures, the solid curves represent the fixed target results, while the dashed curves represent the moving target results. We see that the inertial coupling affects the returning probability only slightly on the subpicosecond time scale. Figure 5 displays the inertial coupling effect on the rate coefficient for each of the four cases considered in Figure 4. In all cases, we can hardly notice the deviation of the fixed target result (solid curve) from the moving target result (dashed curve).
We generalized the BD method recently proposed23 for calculating the time-dependent rate coefficient of diffusioninfluenced bimolecular reactions to include the inertial dynamic effect. Although a description has been given for the case involving simple spherical reactant molecules with isotropic reactivity, the method can be easily generalized to reactions of a more complicated type involving reactants of arbitrary shape and interaction potential. However, usability of the present BD method is restricted to the low reactant concentration case. The reactive superposition approximation given by eqs 2.9 and 2.10 neglects the like-particle correlation effects such as the reactioninduced segregation of dissimilar reactant molecules that may develop at very long times and the excluded volumes among like reactants that can be important when the concentration is high. The present BD method enables us to estimate the rate coefficient very efficiently by using only about 10 000 trajectories. It requires less computing time than numerically solving the Fokker-Planck equation for the relative motion of reactants by using the finite element method. Another key advantage of the present BD method is that it requires only one set of BD simulations to calculate the rate coefficients for any value of the intrinsic rate constant. The numerical method for directly solving the Fokker-Planck equation, which consumes more computing time even for a single calculation, requires separate calculations for different values of the intrinsic rate constant. In the present work, we investigated the inertial dynamic effect on the diffusion-influenced reaction rate by using the proposed BD method. For a reaction between two identical hardsphere reactants, we found that the reaction rate decreases as the mass of the reactant increases but with the product of its size and diffusion constant fixed. When the reactants do not have identical masses and friction coefficients, it was indicated22 that there could be a subtle inertial dynamic effect that arises from the coupling between the CM motion and the relative motion of the reactants. However, we found that the inertial coupling effect should be practically negligible. Acknowledgment. This work was supported by Grant 30520000005 from the Basic Research Program of the Korea Science & Engineering Foundation. References and Notes (1) Smoluchowski, M. V. Ann. Phys. 1915, 48, 1103. (b) Z. Phys. Chem. (Leipzig) 1917, 92, 129. (2) Calef, D. F.; Deutch, J. M. Ann. ReV. Phys. Chem. 1983, 34, 493. (3) Goesele, U. M. Prog. React. Kinet. 1984, 13, 63. (4) Rice, S. A. ComprehensiVe Chemical Kinetics 25. Diffusion-Limited Reactions; Elsevier: Amsterdam, 1985. (5) Ovchinnikov, A. A.; Timashev, S. F.; Belyy, A. A. Kinetics of Diffusion Controlled Chemical Processes; Nova Science: Commack, Russia, 1989. (6) Kotomin, E.; Kuzovkov, V. ComprehensiVe Chemical Kinetics 34. Modern Aspects of Diffusion-Controlled Reactions; Elsevier: Amsterdam, 1996. (7) Montroll, E. W. J. Chem. Phys. 1946, 14, 202. (8) Friedman, H. L. J. Phys. Chem. 1966, 70, 3931. (9) Solc, K.; Stockmayer, W. H. J. Chem. Phys. 1971, 54, 2981. (b) Int. J. Chem. Kinet. 1973, 5, 733. (10) McCammon, J. A.; Northrup, S. H. Nature (London) 1981, 293, 316. (11) Dong, W.; Baros, F.; Andre, J. C. J. Chem. Phys. 1989, 91, 4643. (12) Dong, W.; Andre, J. C. J. Chem. Phys. 1994, 101, 299.
6024 J. Phys. Chem. B, Vol. 105, No. 25, 2001 (13) Zhou, H.-X.; Szabo, A. J. Chem. Phys. 1991, 95, 5948. (14) Jung, Y.; Lee, S. J. Phys. Chem. A 1997, 101, 5255. (b) Lee, J.; Sung, J.; Lee, S. J. Chem. Phys. 2000, 113, 8686. (15) Ibuki, K.; Ueno, M. J. Chem. Phys. 1997, 106, 10113. (b) Bull. Chem. Soc. Jpn. 1997, 70, 543. (16) Harris, S. J. Chem. Phys. 1981, 75, 3103. (b) 1982, 77, 934. (c) 1983, 78, 4698. (17) Naqvi, K. R.; Mork, K. J.; Waldenstrøm, S. Phys. ReV. Lett. 1982, 49, 304. (b) Naqvi, K. R.; Waldenstrøm, S.; Mork, K. J. J. Chem. Phys. 1983, 78, 2710. (18) Sceats, M. G. J. Chem. Phys. 1986, 84, 5206. (b) J. Colloid Interface Sci. 1989, 129, 105. (19) Molski, A. Chem. Phys. Lett. 1988, 148, 562. (20) Chandrasekhar, S. ReV. Mod. Phys. 1943, 15, 1. (b) Wang M. C.; Uhlenbeck, G. E. ReV. Mod. Phys. 1945, 17, 323.
Yang et al. (21) Risken, H. The Fokker-Planck Equation; Springer-Verlag: Berlin, 1989. (22) Rodger, P. M.; Sceats, M. G. J. Chem. Phys. 1985, 83, 3358. (23) Yang, S.; Kim, J.; Lee, S. J. Chem. Phys. 1999, 111, 10119. (24) Lee, S.; Karplus, M. J. Chem. Phys. 1987, 86, 1883. (b) 1987, 86, 1904. (25) Monchick, L.; Magee, J. L.; Samuel, A. H. J. Chem. Phys. 1957, 26, 935. (b) Waite, T. R. Phys. ReV. 1957, 107, 463. (26) Ermak, D. L.; Buckholz, H. J. Comput. Phys. 1980, 35, 169. (27) IMSL Library Reference Manual, Ver. 1.1; IMSL: Houston, 1989. (28) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (29) Szabo, A. J. Phys. Chem. 1989, 93, 6929.
