6430
J . Phys. Chem. 1991, 95,6430-6437
ARTICLES Molecular Dynamics Study of Dlpolar Relaxation P. VijayaKumar School of Chemistry, University of Hyderabad, Central University P.O., Hyderabad 500134, India
and B. L. Tembe* Department of Chemistry, Indian Institute of Technology, Powai, Bombay 400076, India (Received: October 30, 1990: In Final Form: March 20, 1991)
The quilibrium and nonequilibrium solvation dynamics in a nonspherical and nonpolarizable dipolar liquid around charged solute particles (separated by a distance r = 5 A) is studied by using molecular dynamics. The memory kernel for dielectric friction of the solvent obtained from forceforce autocorrelations is analyzed within the framework of Grote-Hynes theory for activated barrier crossing kinetics. The relaxation times of the polarization fluctuations at and around the reactant sites which are needed in the theory of electron transfer proposed by Tembe et al. are evaluated. Solvent dynamics during the process of dielectric relaxation is also studied by using the microscopic and macroscopic dipoledipole time correlation functions for both equilibrium and nonequilibrium solvation around the reactant ion pair. The nonequilibrium dynamics is also discussed in terms of time-dependent polarization fluctuations in different regions of solvent.
I. Introduction Structure and dynamics in dipolar liquids have been subjects of great interest in theoreticaP6 as well as e~perimentalI~-~O investigations. The details of the dynamics of the time-dependent solvent response to a sudden redistribution of charges in a polar solute molecule or among charged reactants have become essential for understanding the role of solvent in studies of charge-transfer processes. In most of the experimental and theoretical studies, attempts were made to create a situation where the dynamics of solvation can be studied in the absence of any chemical charge transfer. Recently, some of the time-resolved studies of solvation of newly created ions and dipoles in dipolar media have been (1) Nee, T. W.; Zwanzig, R. J . Chem. Phys. 1970,52,6353. (2) (a) Wolynes, P. G. J. Chem. Phys. 1978,68,473; (b) A m . Rm.Phys. Chem. 1980,31,345; (c) J . Chem. Phys. 1987,86, 5133. (3) (a) Calef, D. F.; Wolynes, P. G. J . Chem. Phys. 1983, 78, 4145; (b) J . Phys. Chem. 1983,87, 3387. (4) Bagchi, B.; Oxtoby, D. W.; Fleming, G. R. Chem. Phys. 1984,86,257. (5) van der Zwan, G.; Hynes, J. T. J . Phys. Chem. 1985, 89, 4181. (6) Loring. R. F.; Yan, Y. J.; Mukamel, S. J . Chem. Phys. 1987,87,5840. (7) Bagchi, B.; Castner, E. W., Jr.; Fleming, G. R. J . Mol. Srrucr. 1989, 194, 171.
(8) Castner, E. W., Jr.; Fleming, G. R.; Bagchi, B.; Maroncelli, M. J . Chem. Phys. 1988,89, 3519. (9) Nichols 111, A. L.; Calef, D. F. J . Chem. Phys. 1988,89, 3783. (IO) Bagchi. B.; Chandra, A. Proc. Indiun Acud. Sci. Chem. Sci. 1988, 100,353; J . Chem. Phys. 1989.90,7338; Chem. Phys. I r t r . 1989,155,533. (1 1) Chandra, A.; Bagchi, B. Chem. Phys. Lcrr. 1988, 151,47; J . Phys. Chem. 1989, 93, 6996. (12) Bagchi, B. Annu. Reu. Phys. Chem. 1989, 40, 115. (13) Carter, E. A.; Hynes, J. T. J . Phys. Chem. 1989, 93, 2184. (14) Zichi, D. A.; Ciccotti, 0.;Hynes, J. T.; Ferrario, M. J . Phys. Chem. 1989, 93, 6261. (15) Bagchi, B.; Chandra, A.; Fleming, G. R. J . Phys. Chem. 1990, 94, 5197. (16) Wei, D.; Patey, G. N. J . Chem. Phys. 1990.93, 1339. (17) Maroncelli, M.; Fleming, G. R. J . Chem. Phys. 1987, 86. 6221. (18) Simon, J . D.;Su, S.-G. J . Chem. Phys. 1987,87, 7016. (19) Su, S.-G.; Simon. J. D. J . Chem. Phys. 1988, 89, 908. (20) Jarzeba, W.; Walker, G. C.; Johnson, A. E.; Kahlow, M. A.; Barbara, P. F. J . Phys. Chem. 1988, 92, 7039.
0022-3654/91/2095-6430$02.50/0
rep~rted.'~-~O In another investigation, theoretical and experimental results of solvation dynamics around a newly created charge were compared.21v22 In an analysis of solvent response, it was first suggested by Boyd that there will be a lag in the ability of the solvent to readjust to the new environment due to the dielectric friction experienced by the solvent molecules. The dielectric Boyd,u Zwanzig,25 friction picture was extensively used by FUOSS,~ and Hubbard and Onsager.26 In almost all the theories the total friction coefficient ttohl is separated into a viscous part, tv, and an electrical contribution or dielectric friction, tD:
= t V + tD (1) where tvis given by Stokes law: tv = 4uqa (2) Here 7 is the viscosity, and a the radius of the particle. The time-dependent dielectric friction fD(t) plays a very crucial role in the theory of chemical reaction rates in condensed phases. Grote and Hynes (GH)proposed a rate equation for the activated barrier crossing kineticsz7by extending Kramers' theory with the introduction of time-dependent friction (TDF). In Kramers' theory?* t ( t ) is assumed to be equal to r(t) = to6(t), where 6 ( t ) is the Dirac 6 function and tois the friction coefficient. In G H theory, it is further assumed that the memory kernel is not influenced by the activation barrier for the reaction, which means that one can evaluate TDF in the absence of any barrier and still apply the same function to processes where a barrier is involved. ftotal
(21) Karim, 0. A.; Haymet, A. D. J.; Banet, M. J.; Simon, J. D. J . Phys. Chem. 1988, 92, 3391. (22) Maroncelli, M.; Fleming, G. R. J . Chem. Phys. 1988,89, 875. (23) Fuoss, R. M. Proc. Narl. Acud. Sci. W.S.A. 1959, 45, 807. (24) Boyd, R. H. J . Chem. Phys. 1961.35, 1281; J . Chem. Phys. 1963,
-39.- (25) .--2376.-(a)- Zwanzig, R. W. J . Chem. Phys. 1961, 35, 1281; (b) J . Chem. Phys. 1963,38, 1603; (c) J . Chem. Phys. 1%3,38,1605; (d) J . Chem. Phys. 1970. -.. -, 52. - -, 3625. - - -- . (26) Hubbard, J.; Onsager, L. J . Chem. Phys. 1977,67, 4850. (27) Grote, R. F.; Hynes, J. T. J . Chem. Phys. 1980, 73,2715; 1981, 75, 2191. (28) Kramers, H. A. Physica 1940, 7 , 284.
0 1991 American Chemical Society
The Journal of Physical Chemistry, Vol. 95, No. 17, 1991 6431
Molecular Dynamics Study of Dipolar Relaxation The dynamical behavior of the solvent can be monitored by calculating the time-dependent cavity field (the electric field in the neighborhood Pf the reactants) resulting from the polarization fluctuations in the medium, The time-dependent solvent polarization is related to the time=dependent solvation energy whose effects are observed in the time-dependent fluorescence experiments. In most of the experimental studies, the polarization fluctuations are studied by using the normalized spectral shift time correlation fuqction. In theoretical studies, the time-dependent where ?L is polarization, P(r,t), is dccorflposed as 9 = 3 L + the longitudinal part and PT is the transverse part.29-30 The relaxation of the transveq part is the same as the relaxation qf the displacement vector, D,after a jump in the electric field, E. This process is associated wit! the Debye relaxatgn time TD. The relaxation of electric field, E, after a jump in D is a relaxation of longitudinal part of the polarization vector field, and this is associated with the longitudinal relaxation time 7L. The polarization fluctuations in our simulations are studied by observing the time-dependent cavity field fluctuations at and around the solute ions. The outline of the paper is as follows. In section I1 we discuss the theories of dielectric friction coefficient and cavity field fluctuations used in our calculations. We also discuss the use of TDF in the activated barrier crossing charge-transfer processes followed by the relevant cavity field time correlation functions (CFTCFs) in a theory of electron-transfer reactions. The details of the model and the method used in our simulations are described in section 111. In section IV we present the results and discussions of the TDF memory kernels and compare these with continuum theory results. Then we present our results of the CFTCFs, followed by a description of the dielectric relaxation with the help of the macroscopic and microscopic dipole moment correlation functions. In the last part of section IV, we present our results of nonequilibrium dynamical solvation around the newly created charge distribution. The summary and conclusions of the present study are given in section V.
9,
11. Theory
The dielectric friction coefficients lD’s for an ion and a dipole in a continuum theory are given by the expressions25w (3)
m dr(t)/dt = F(f) -
-Sy
-
-x‘{(t
- s) V(s) ds
(5)
Grote and Hynes introduced a frequency-dependent friction, {(a)into the Kramers’ theory of activated barrier crossing in
condensed phases to account for the non-Markovian effects. In their theory the dynamics is described by the generalized Langevin equation2’ (29) Kivelson, D.; Friedman, H. L. J . fhys. Chem. 1989, 93, 7026. (30) Newton, M. D.; Friedman, H. L. J . Chem. fhys. 1988, 88, 4460. (31) Nakahara, M.; Ibuki. K. J . fhys. Chem. 1986, 90, 3026.
(6)
S)
0
where m is the reduced mass, F ( f ) = -aU/dr, R(r) is the random force, and f(s) is the TDF coefficient. The TDF coefficient is related to the temperature and the random force autocorrelation function by the second fluctuation dissipation theorem. In the case of a large particle the random force can be replaced by the systematic ion-solvent force. Then the friction coefficient takes the following { = (1/3knJmdf
(R(t)*R(O))
-
{=
(1/3kr)J-dr
(F(t)-F(O)) (7)
The memory kernel or the time-dependent friction, { ( t h is l(t) = (1/3kT)(F(W(O))
(8)
A full knowledge of the memory kernel is necessary to describe chemical reactions involving a crossing over an activation barrier. The TDF coefficient { ( t ) obtained from eq 8 contains information about both translational (diffusive) and rotational drags, due to dielectric as well as Stokes frictions. The torque on a rotating dipole in a cavity may be expressed in terms of the rotational friction coefficient Crot(t) in a manner analogous to eq 5 by replacing the translational velocity by angular velocity:’
N(t) =
-it{& - s) ds + random torque Q(s)
(9)
In the absence of translational diffusion, the rotational part of memory kernel {rot(t) may be obtained from the torque-torque autocorrelation function under the conditions when eq 8 is valid: {rot(t)
= (1/3kT)(N(t).N(O))
(10)
The TDF has been found to be very useful in evaluating the rate constants for the activated barrier-crossing p r o c e ~ s e s . ~ ~ - ~ ~ One of the widely used rate theories for activated barrier crossing is the Grotc-Hynes (GH) theory.27 They have used TDF in the equation for the rate constant for the first time and modified Kramers’ theory for activated barrier crossing. In their formulation, the rate equation is given by
k = kmT(Xr/@b)
(4)
where to and c, are the static and high-frequency dielectric constants, respectively, q is the charge of the ion, and p is the dipole moment. The importance of time-independent dielectric friction in the theory of ionic mobility was discussed by Hubbard and Onsager.% In an extension of the work of Hubbard and Onsager, formulas were obtained for the residual friction, tD = - tV by fitting the numerical solution of the appropriate electrohydrodynamic equation of m ~ t i o n . ~ ’ The friction term appears in the Langevin equation as a contribution to the retarding force, -Sy, where V is the velocity of the Brownian particle. When there is a lag between the motion of the Brownian particle and the response of the solvent or the heat bath, then the frictional force can be generalized to include a time-dependent friction or the memory kernel { ( t ) :
sh, S(s)v(t - + R(t)
(11)
where
here kTSTis the transition state theory rate constant, A, is the reactive frequency, a b is the barrier frequency, and f(X,) is the Laplace transform of the TDF at frequency A,. In an alternative theory of electron-transfer reaction, the solvent dynamical effects are discussed in terms of cavity field ( G ( t ) ) fluctuations. The influence of the solvation dynamics on the rates of the charge-transfer processes can be estimated from the following relations (taken from the theory proposed by Tembe et al.” and Friedman and Newton36) kAB
xm4*9gAB(r)
R.4B(r>
dr
(12)
(32) Berkowitz, M.; Wan, W. J . Chem. fhys. 1987, 86, 376. (33) Zhu, S.-B.; Lee, J.; Robinson, G.W. J. Chem. Phys. 1988,88,7088. (34) Zhu, S.-B.; Lee,J.; Robinson, G. W.; Lin, S.H. J. Chem. fhys. 1989, 90, 6340. (35) Samanta, A.; Ghosh, S. K.; Sadhukhan, H. K. Chcm. fhys. k r r . 1990. 168,410. (36) Friedman, H. L.; Newton, M. D. Faraday Discuss. Chem. Soc. 1982, 74, 73. (37) Tembe, B. L.; Friedman, H. L.; Newton, M. D. J . Chem. fhys. 1982, 76, 1490.
6132 The Journal of Physical Chemistry. Val. 95, No. 17, 19'91
where
VijayaKumar and Tcmbe = -0Se and +0.5e, respectively). The total solute-solvent interaction potential is given by r
Here k , , ~is the bimolecular rate constant between Fe2+and Fe3+ is the equilibrium pair correelectron-exchangc reaction, gAB(r) lation function, k,,g(r) is the local rate canstant at the separation r, 19 = l/kBT, E%) is the activation energy, 7a is the chariyteristic time for the electron to jump from Fe2+to Fe3+when C ( t ) = 0, and is the relaxation time for the CFTCF. We have studied the behavior of the time correlation functions for the fluctuating electric fields at and in the vicinity of the two reactant cavities. The field at the reactants is closely related to the cavity field of the continuum theories, except that in our simulations the cavity consistsof the region containing two reactant sites. It is shown that if we approximate the medium to be a continuum and assume a Debye form for the dielectric function, Le., e(w) = eo [(b- e.,)/( 1 i w D ) ] ,then the time correlation function for the cavity field takes the following expression?* which is also a single exponential:
+
+
where
Here R is the radius of the cavity. This cavity field is a measure of polarization fluctuations of the medium. The cavity field relaxation time is equivalent to the longitudinal relaxation time f L , if the cavity containing the point dipole is assumed to have a dielectric constant equal to e,. This f 0 may be used as a measure for monitoring the solvent relaxation process. The correlation time of the cavity field will strongly influence the rate of a chargetransfer process in the cavity. We also discuss the dielectric relaxation in terms of microscopic and macroscopic dipole-dipole correlation functions. The molecular rotational relaxation is studied by using the single-molecule dipole-dipole autocorrelation function:
440 =
(cC'(t).cC'(O))
/ (r'(OY)
(16)
Knowledge of the total dipole moment autocorrelation function a(?),which gives information about the collective orientational ordering of the bulk medium, was also considered in the present study:
where M(t) = C f p f ( t ) . 111. Model and Method The system used for simulations consists of 98 dipoles confined to a cubic cell of constant volume with periodic boundary conditions. The solvent dipoles consist of two Lennard-Jones (LJ) centers separated by 1.5 A, and the two LJ centers in each dipole carry charges of -0.25e and +0.25e. This model dipole corresponds to a dipole moment of 1.8 D. The packing fraction of the simulation system is 0.41, which corresponds approximately to a dense dipolar medium. There are no intramolecular interactions, and the intermolecular interactions are modeled by site-site Lennard-Jones interactions and Coulombic interactions:
4f)(ra,3)= #LJ(ra,3) + &(ra& (18) where a and /3 are the sites on molecules i and j , respectively. The solute particles are modeled as LJ spheres separated by 5 A and are fixed in the middle of the simulation box along the x axis. The two solute particles carry an equal and opposite charge (q (38) Hubbard, J. B.;Wolynes, P. G. J . Chem. Phys. 1978, 69, 998.
where qf and qj refer to the charges on the solute ion and the solvent site, respectively, and ru is the distance betwaen the salute ion and the solvent site. All the interactions are treated by the minimum image scheme.3g This scheme is considered so as to increase the number of intermolecular interactions that are explicitly included in this method. In a recent study by Huston and RosskyM it was reported that the minimum image prescription gives satisfactory results, and they also reported that spherical truncation gave some unphysical results in their study. Zhu and Robinson4' have already reported that the absence of solutesolute interactions assumes very dilute concentrations. The fact that the change in system size from 98 to 198 dipoles had no significant effect on the simulation results further supports this assumption. Even in the studies of ionic conductivities of Na+ and C1- ions in aqueous solutions, the minimum image method has been used in estimating condcutivities at infinite dilution.32 The separation between the solutes is kept constant throughout the simulation. The solute is modeled in this way so that it represents a very likey configuration for an electron-transfer reaction in the dipolar medium. We have used the parameters of argon for our LJ centers on solvent molecules as well as for solute particles. All the simulations were performed at 298 K. The dynamics was carried out with the use of the Verlet algorithm42and the SHAKE^^ procedure. The time step used in our simulations is 2 X s. After an equilibration of 15K time steps, we have performed the simulation for 20K time steps to study the equilibrium dynamics. We evaluated the electric fields at and around the solute ions and the single molecule (microscopic) and total system (macroscopic) dipole moments during the simulation after every time step. We also store the forces acting on the solute ions after each time step in the simulation and also the configurations after every 2 ps. Each one of these configurations provides the input for the simulations in which we studied the nonequilibrium dynamics. We examined the nonequilibrium dynamics by introducing a perturbation in the equilibrated system brought about by interchanging the charges on the solutes. We use this perturbed configuration and carry out the simulation for 4 ps starting from each one of the above 20 configurations. This enables us to study the nonequilibrium dynamics. The model system that we investigated here can be used to examine the validity of various continuum and molecular theories for dielectric relaxation and solvation dynamics. To check the finite size effects, we performed a test simulation with 198 dipolar molecules at the same density and compared these results with the simulation employing 98 dipoles. The differences between these two results are within 5%. This indicates that the increase in the size did not effect the dynamics of the system. Similar observations have been made in solvent relaxation studies by Engstrom et al.," who studied systems with 64 and 125 water molecules, and also by Karim et al.,2' who studied systems with 212 and 504 water molecules. In another test, we performed a simulation with 98 dipoles where we made use of the reaction field (RF) approximation to incorporate the polarization of the continuum beyond the cutoff distance R,, induced by the net dipole inside the sphere of radius R, centered on that molecule.45 The reaction field adds a term to the potential energy of a charged (39) Allen. M. P.; Tildesley, D. J. Computer Sirnufatiom of Liquids; Clarendon Press: Oxford, 1987. (40) Huston, S. E.; Rcwsky, P. J. J . Phys. Chem. 1989, 93, 7888. (41) Zhu. S.-B.; Robinson, G. W. J . Chem. Phys. 1989, 90. 7127. (42) Verlet, L. Phys. Rev. 1967, 159, 98. (43) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J . Compur. Phys. 1977, 23, 327. (44) Engstram, S.;J h n . B.; Impey, R. W. J . Chem. Phys. 1984, 80, 5481. (45) Bottcher, C. J . F.; Bordewijk, P. Theory of Electric Polarization; Elsevier: Amsterdam, 1978.
Molecular Dynamics Study of Dipolar Relaxation
Figure 1. Normalized time-dependent friction coefficients calculated from force-force autocorrelation functions at the -ve ion (a) and at the +ve ion (b).
site i due to all other charged sites j within the cutoff sphere R, centered about the site i.& The additional term in the potential is given by
In our simulations with reaction field method, conducting boundary conditions were used with t R F = m, and the cutoff distance is half the length of the box. We observe about 10% differences in the results obtained in this case. This difference is essentially due to the additional term in the potential given by vRF(i) of eq 20. The study by Neumann et ala4' indicates that the reaction field method and the method of Ewald sums for treating the boundary effects yield similar results in the estimation of dielectric properties. In another study by Gray et it was shown that for N L 108 (Nis the number of molecules in the simulation cell), both reaction field and Ewald method have converged to the infinitesystem limit. We conclude that the minor differences in the results obtained for a larger system and the system with reaction field do not affect the conclusions of the present study.
IV. Results and Mscussion (A) Time-Dependent Friction. The dynamical characteristics of the motion of the particles surrounding the solvated ions can be characterized by the TDF coefficient, { ( t ) . A microscopic expression for this drag coefficient is given in terms of the correlation function of the random forces exerted on the solute molecules. In the case of a large solute ion, this is the same as the autocorrelation function of the force exerted on it due to regular ion-solvent interactions (eq 7). This approximation can also be extended further to the case of a fured ion to evaluate K t ) from the regular ionsolvent force-fm autocorrelation function.2. We numerically evaluate the memory kernel for the friction coefficient at the two solute (reactant) cavities in our simulations. Our memory kernels include correlations of both the hard (H) and the soft (S)parts of the forces and the force-force autocorrelation function has the following form: (F(t).F(O)) = [(FH(t)'FH(t)) + (FH(t).Fs(t)) + (Fs(r)*FH(t)) + (Fs(Ws(t))1 (21) The memory kernel curves calculated in our simulations are shown in Figure 1 . It is seen from these curves that these (46) Ncumann, M. Mol. Phys. 1983,50, 841. (47) Ncumann, M.; Steinhauser. 0.;Pawley, G. S. Mol. Phys. IW, 52,
The Journal of Physical Chemistry, Vol. 95, No. 17, 1991 6433
Figure 2. Normalized time-dependent friction coefficient f i t ) calculated from the continuum model for water (a), methanol (b), model I (c), model I1 (d), and model I11 (e). See Table I for the parameters used in obtaining these curves.
functions deviate from the monoexponential behavior. The TDFs at the two solute ions are identical up to 0.5 ps, and the differences between the two beyond 0.5 ps are within the percentage errors in the simulation. This is expected in our system as long as the charges on the two solute sites are equal in magnitude. In our simulations we have a solute ion pair A0.S--.B0.5+separated by a fixed distance of 5 A, which may be assumed to be a representative of a dipolar molecule inside a cavity surrounded by the solvent. The charge-transfer process corresponds to a flip of the dipole from A0.5--.B0.5+to A0.S+-.B0.5-.This is equivalent to a rotation of the dipole inside the cavity. The TDF coefficients obtained from the equilibrium simulations in the present work give information about the solvent response that would result from small variations in the solute ion's charge distribution as long as the variations are in the linear response regime. This enables us to compare our TDF coefficients with the Kt)'s obtained for a model for dipole i s o m e r i z a t i ~ nwhich ~ ~ ~ ~are shown in Figure 2. This comparison gives us some indications as to which actual molecular liquid our model dipolar medium corresponds. The expression for {(t), i.e., the time-dependent rotational dielectric friction coefficient for the dipole (in the presence of spatial diffusion), is given by
where r
1 €0
-1 to + €,WTD 2 1 +UTD
+ D,R + (1/3)D?R2 + 1 + D.R 1 + D1R + (2/3)DI2R2+ (1/3)D13R3
11
1 + DIR
and D, and Dl are defined by D, = [ ( t o D ~ ~ ) - ' (+c o€ , 0 7 ~ ) ] ~ / *
1
DI = [(%)-'(I + W T D ) ] ' / * here T L = [(2t, + I ) / ( & + l)]rD, is the longitudinal relaxation
_..
07
(48) Gray, C. 0.;Saingcr,Y.s.; Jcelin. C. G.; Cummings, P.T.; Goldman, S. J. Chem. Phys. 1986,85, 1502.
(49) van der Zwan, G.; Hynes, J. T. Physica 1983, IZIA, 228. (50) van der Zwan, G.; Hynes, J. T. Chem. Phys. &If. 1983, 101, 367.
6434 The Journal of Physical Chemistry, Vol. 95, No. 17, 1991
VijayaKumar and Tembe
TABLE I: Dielectric Properties,Relaxation Times, and Diffusion Coerficients~
solvent
CHiOH H20 model I
model 11 model 111
CO
6-
TD, ps
TL, ps
D,A2/ps
32.62 78.6 78.6 32.62 78.6
5.6 1.85 1.85 5.6 1.85
47.7 8.6 8.6 47.7 8.6
8.22 0.24 0.24 8.22 0.3
0.294 0.25 0.39 0.39 0.39
OSee ref 5 for dielectric data and relaxation times. time and L-‘represents the inverse Laplace transform. The friction at r = 0, Le., {(PO) is
where I is the moment of inertia, R is the radius of the spherical cavity, fi is the dipole moment, D is the diffusion coefficient, and 7D is the Debye relaxation time. In Figure 2 we have shown the TDF curves generated within the continuum theory framework for water (curve a) and methanol (curve b) and also the TDF plots (curves c-e) generated by using the parameters of our dipolar model. The parameters used in evaluating f(t) of Figure 2 are shown in Table I. For curves c and d, the values of molecular parameters cc and I are the same as those in our model potential (eq 18), and the value of the diffusion coefficients is from the plots of mean-square displacements in our simulations. Since we have not estimated the dielectric constants of our model, we have used the values of water and methanol to obtain curves c and d. It is observed from the graphs in Figure 2 that curves a and c relax on a shorter time scale than the curves b and d. The reason for this can be attributed to the values used for the parameters 6.c,, 7L,and 7D in evaluating eq 22. It is known that for liquids with large values of e,, and small values of c, the decay of { ( t ) is rapid? We also find that curve c relaxes faster than curve a. This is due to the larger value (Le., the value obtained in the simulations) for the diffusion coefficient used in obtaining curve c than the value used to evaluate curve a. The same reason accounts for the faster decay of curve d than the decay of curve b. We find that, although the initial parts of f(t) evaluated using continuum theory (curves a and c) resemble the s‘(f)’s obtained from our simulations (Figure l), they are found to decay at a slightly faster time scale. This indicates that the 7L value we used in evaluating curve c is rather small. Indeed, if we use the value of T~ = 0.3 ps, we obtain curve e. This value of 0.3 ps was obtained from the C F E F s , which will be discussed later in this section. We find that curve e (Figure 2) resembles closely the memory kernels evaluated from our simulations shown in Figure 1. Although the { ( t ) ’ s obtained from our simulations (Figure 1) resemble waterlike models at short times, at longer times there is a deviation in their behavior. This may be due to the molecular nature of the solvent in our simulations. This behavior may also be due to the larger size of the model solvent that we used in our simulations. It is observed that as the molecular size increases, these curves decay more slowly. The similarity in the behavior of TDF of our model (Figure 1 ) and the TDF of water indicates that the dielectric properties (Le., e’s, 7L and T ~ of) our model resemble the dielectric properties of water. This is further corroborated by the fact that the relaxation rates of the TDF curves depend on the static and optical dielectric constants as weK9 The exact behavior of the TDF as a function of these dielectric constants and the shape of the molecule needs to be investigated in greater detail. (B) Cavity Field Fluctuatiom. To study the solvation dynamics due to polarization fluctuations of the medium in our dipolar system, we evaluated the time-dependent cavity fields at the solute sites and at the center of the system due to the solvent after each dynamics time step. We examine the following time correlation function to investigate the details of the dynamics around the solute ions: C(f) = (&f).G(O))/(G(0)2)
(24)
‘
-0.4
1
I
O h 0 025 0.50
1
0.75 1.00 1.25 1.50 1.75 2 10 Time (ps)
Figure 3. Cavity field time correlation functions at the -ve ion (a), at the +ve ion (b), and at the center of the system (c).
where 6(t)is the time-dependent total electric field at the ions or at the center of the simulation cell. The electric fields a t the reactant sites are the same as the cavity fields. In our simulations we can think of our reacting system as a nearly ellipsoidal cavity of low dielectric constant tCimmersed in a frequency-dependent dielectric medium. The dynamical effects that correspond to Onsager’s reaction field in a continuum model were discussed by Nee and Zwanzig.’ There, a rotating dipole in the cavity experiences a lagging time-dependent reaction field from the medium that tends to slow down the rotating dipole. Hubbard and Wolynes introduced the concept of the dynamLcal cavity field.38 We have calculated the dynamical cavity field G(t) in our simulations. The plots of CFTCFs at the -ve and the +ve ions and a t the center of the system are shown in Figure 3. The cavity field is a measure of the polarization fluctuations of the medium that may be estimated in linear response theory in a medium characterized by . relaxation a frequency-dependent dielectric function ~ ( u )The times ( 7 6 ) for the CFTCFs can be obtained by fitting the initial parts of these curves to single exponentials. However, none of these curves behave monoexponentially at longer times. These T ~ ’ have S an important role in the studies of solvation dynamics and their influence on the rates of the charge-transfer processes. We have fitted the initial 0.25-ps parts of the CFTCFs shown in Figure 3 to single exponentials and found that the relaxation times were falling between 0.25 and 0.3 ps. Continuum theory calculations lead to the value of 76 = 0.3 ps for water at 25 O C . The value for 76 that has to be used in eqs 12 qnd 13 should be obtained from the time integral of the CFTCF:
If CFTCFs attain a plateau value of zero beyond a time 7p, then the upper limit of the integration in the above equation may be taken to be 7p.32In our studies, the plateaus do not seem to have been attained even at T~ = 2 ps; however, an integration up to T~ = 2 ps may not give an accurate estimate of 70 due to the relatively larger errors in CFTCFs a t 7p 1 1 ps. Thus we have used the initial part of the CFTCFs for obtaining an estimate of T ~ Improved . estimates may be obtained by performing a longtime simulation from which it will be possible to obtain the CFTCFs accurately at the longer times. We are investigating a more realistic model for water to make an estimation of the cavity field relaxation times involved in aqueous charge-transfer processes. This will lead to more accurate values for 70. These T~ values can be used in eqs 12 and 13 to see the affect of polarization fluctuations of the medium on the rates of the charge-transfer processes. (C) Dipole-Dipole Correlation Functions. Another convenient approach to the study of polarization relaxations is the calculation
The Journal of Physicql Chemistry, Vol. 95, No. 17, 1991 6435
Molecular Dynamics Study of Dipolar Relaxation
00 ! a0
1
I
I
0
10 1.5 20 Timo(ps1 Figure 4. Normalized autocorrelation functions of the single-molecule dipole moment (a) and total dipole moment of the system (b). 0.3
of the microscopic and macroscopic dipole-dipole autocorrelation functions. It is possible to obtain the frequency-dependent dielectric constant from the total dipole moment autocorrelation function @(r).S1 This correlation function (eq 17) for our system is shown in Figure 4 (curve b). The relaxation time corresponding to this function represents the dielectric relaxation time. This is the same as the time required for the collective orientational reordering of the system. In our simulations we found that the total dipole moment correlation ((M(t).M(O))) function is relaxing faster than those reported in the works of Anderson et al.51and K u ~ a l i k .In ~ ~a recent study of the total dipole moment correlation function which uses a reaction field dielectric constant cRF for estimating the long-range electrostatic forces, it was found that a decrease in the value of this t R F leads to a faster r e l a x a t i ~ n . ~ ~ In our simulations, the same reason can be attributed to the faster relaxation of the @(r), since the minimum image method employed in our simulations uses a value of t = 1 while calculating the interactions in the medium. Another related quantity is the single molecule correlation function 4 ( t ) defined by eq 16. The time TS is the single-molecule correlation time and is defined by Ts
= Jmdt (rl(t).ccl(O))/(r1(0)2)
(26)
Both TS and f D are similar in magnitude, but T~ is smaller than TD. In general, for liquids with large to we have29
