J. Phys. Chem. 1992, 96, 4068-4074
4068
a result that will usually not conform to the experimental evidence. Adding a linear as well as a logarithmic term to g g = s + h/t bT + c In T (~5)
+
we have
AcP = -~1(~2(26T + C)
('46) Thus we see that a linear dependence of Acp on T directly leads
to expression A5 for g. In practice such a Ac,( 7") behavior will represent a first approximation only, but it serves the present qualitative purpose. The g(T) function IOd,used in the text for illustration, thus involves b = 0 and makes Acp nonzero but independent of T . Registry NO. PMMA, 9011-14-7.
Dielectric Friction and Solvation Dynamics: A Molecular Dynamics Study Margaret Bruehlt and James T. Hynes* Department of Chemistry and Biochemistry, University of Colorado, Boulder, Colorado 80309-021 5 (Received: October 17, 1991)
The van der Zwan-Hynes relation connecting the solvation time and the dielectric friction for a solute in a polar solvent is tested via molecular dynamics computer simulation. For partially and fully ionic diatomic solute pairs in a model polar aprotic solvent, for which there is considerable dielectric friction, the relation is found to be satisfied to within a factor of 2 and to correctly follow the trends observed for different solute pairs. Deficienciesof the relation are also discussed. Pronounced solute rotational caging associated with strong electrostatic solutesolvent interactions is also observed.
1. Introduction
The concept of the dielectric friction' on a translating ion or a rotating dipole immersed in a polar solvent has a long history, beginning with the work of Born in 1920.* For example, the time lag of the surrounding polar solvent electric polarization for a steadily rotating dipole leads,3 in a simple dielectric continuum description, to a dielectric friction constant proportional to the longitudinal dielectric relaxation time4 of the solvent. Considerable effort has been made to experimentally test this and related predictions in the context of molecular reorientation An independent measure of the dynamics of a polar solvent in the presence of dipolar (and/or charged) solutes is the so-called 'solvation time" T ~available , from time-dependent-fluorescence (TDF) measurements.* For example, in a Franck-Condon electronic transition in which there is an excited-state dipole moment jie different from the ground-state value iig, there will be subsequent solvent relaxation from the equilibrium appropriate to ii,to that appropriate to jie; this relaxation can be related to the time evolution of the TDF Stokes shift.8-'0 A convenient molecular measure of the solvation time scale is the time area of the normalized equilibrium time correlation function (tcf)
should approximately hold, connecting the solvation time to the ratio of the dielectric friction constant {d to the initial time value
(1) For reviews see: Wolynes, P. G. Annu. Reu. Phys. Chem. 1980, 31, 345; Madden, P. A.; Kivelson, D. Adv. Chem. Phys. 1984, 56, 479. (2) Born, M. Z . Phys. 1920, 1 , 221. (3) Nee, T.; Zwanzig, R. J . Chem. Phys. 1970,52,6353. Hubbard, J. B.; Wolynes, P. G. J. Chem. Phys. 1978, 69, 998. (4) See,e.g.: Hubbard, J.; Onsager, J. J . Chem. Phys. 1977, 67, 4850; Hubbard, J. B. J . Chem. Phys. 1978, 68, 1649. ( 5 ) (a) Templeton, E. F. G.; Kenney-Wallace, G. A. J. Phys. Chem. 1986, 90, 2896, 5441. (b) Blanchard, G. J. J . Phys. Chem. 1988, 92, 6303. (c) Kivelson, D.; Spears, K. G. J. Phys. Chem. 1985, 89, 1999. (d) Phillips, L. A.; Webb, S.P.; Yeh, S. W.; Clark, J. H. J . Chem. Phys. 1985,89, 17. (e) Alavi, D. S.; Waldeck, D. H. J . Phys. Chem. 1991,95,4848. ( f ) Nakahara, M.; Ibuki, K. J . Chem. Phys. 1986,85,4654. (g) Ben-Amotz, D.; Scott, T. W. J. Chem. Phys. 1987,87, 3739. (6) Simon, J. D.; Thompson, P. A. J . Chem. Phys. 1990, 92,2891. (7) Alavi, D. S.; Hartman, R. S.; Waldeck, D. H. J . Chem. Phys. 1991, 94,4509. (8) For r a n t reviews, see, e.g.: (a) Barbara, P. F.; Jarzeba, W. Ado. Photochem. 1990, 15, 1. (b) Barbara, P. F.; Kang, T. J.; Jarzeba, W.; Fonseca, T. In Perspectives in Photosynthesis; Jortner, J., Pullman, B., Eds.; Kluwer: Dordrecht, 1990; p 273. (c) Simon, J. D. Acc. Chem. Res. 1988, 21, 128. (d) Kosower, E. M.; Huppert, D. Annu. Rev. Phys. Chem. 1986,37, 127. (e) Bagchi, B. Annu. Reu. Phys. Chem. 1989,40, 115. ( f ) Maroncelli, M.; MacInnis, J.; Fleming, G. R. Science 1989, 243, 1674. (g) Maroncelli, M. J . Mol. 139. 1991, in press. (9) Bagchi, B.; Oxtoby, D. W.; Fleming, G. R. Chem. Phys. 1984,86,257. (10) van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 1985,89,4181. Note that the friction defined in this reference is equal to the product of the friction Here 6 A E is the fluctuation in the difference, for fixed solvent and the moment of inertia I in the present work. molecular coordinates, of the solute-solvent interaction energy (11) (a) Carter, E. A.; Hynes, J. T. J . Chem. Phys. 1991, 94, 5961. (b) for the excited and ground dipole The molecular See also: Hynes, J. T.; Carter, E. A.; Ciccotti, G.; Kim, H. J.; Zichi, D. A.; formula eq 1.1 can be identified with the experimental solvation Ferrario, M.; Kapral, R. In Perspecriues in Photosynthesis; Jortner, J., Pullman, B., Eds.; Kluwer: Dordrecht, 1990; p 133. time T~ mentioned above, upon the assumption of linear response: (12) Maroncelli, M. J . Chem. Phys. 1991, 94, 2084. T~ = T ~ There ~ can 'be some ~ deviations from a strict validity (13) Maroncelli, M.; Fleming, G. R. J. Chem. Phys. 1988, 89, 5044. of this a s s u m p t i ~ n , l ' - ~but ~ * 'these ~ will not directly concern us (14) Bader, J. S.; Chandler, D. Chem. Phys. Lett. 1989,157, 501. Levy, here. R. M.; Kitchen, D. B.; Blair, J. T.; Krogh-Jespersen, K. J . Phys. Chem. 1990, 94, 4470. These two solvent dynamical measures have been approximately (15) Fonseca, T.; Ladanyi, B. M. J . Phys. Chem. 1991, 95, 2116. related by van der Zwan and Hynes (vdZ-H),I0 in an effort (16) (a) van der Zwan, G.; Hynes, J. T. J . Chem. Phys. 1982, 76, 2993; designed to provide an experimental route to the independent J . Chem. Phys. 1983, 78, 4174; Chem. Phys. 1984, 90, 21. Zichi, D.A.; prediction of dielectric friction effects in chemical reactions inHynes, J. T. J. Chem. Phys. 1988,88,2513. (b) van der Zwan, G.; Hynes, volving nuclear dipolar realignment or charge di~p1acement.l~~'~ J. T. Chem. Phys. Left. 1983, 101, 367. (c) Ashcroft, J.; Besnard, M.; Aquada, A.; Jonas, J. Chem. Phys. Lett. 1984, 110,420. For the dipolar solute case, these authors showed that the relation (17) (a) Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J. Chem. Phys. 1989, 90, 3537. (b) Bergsma, J. P.; Gertner, B. J.; Wilson, K. R.; Hynes, J. T. J . TS = cd/cd(t=O) (1.2) Chem. Phys. 1987.86, 1356. Keirstead, W.; Wilson, K. R.; Hynes, J. T. J . Chem. Phys. 1991, 95, 5256. Roux, B.; Karplus, M. J. Phys. Chem. 1991, 95, 4856. 'NSF Postdoctoral Fellow 1990-1992.
0022-3654/92/2096-4068$03.00/0
0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 4069
Dielectric Friction and Solvation Dynamics of {&), the time-dependent dielectric friction. The latter is the dielectric component of the friction in the generalized Langevin -.3,10.18 equation (GLE) for the solute angular velocity w. d3( t ) = - i ‘ d T {(t-T) Z ( T ) 4dt
d(t)
(1.3)
where d(t)is a “random” torque, divided by the square root of the solute moment of inertia I. The isolation of this component will be discussed in section 2.A. In the special case of a dielectric continuum treatment in which the solute is treated as a point dipole in a spherical cavity, each side of eq 1.2 would be simply proportional to the solvent longitudinal relaxation time.I0 An important feature of eq 1.2 is, however, that it is not limited to such a continuum approximation. The molecular friction also contains solvent orientational and translational contribution^,^^,'^,^^,^^^ as does the molecular solvation tcf, in eq 1.1. The basic vdZ-H relationlo has recently been tested experimentally by Simon and Thompson6 and by Waldeck and coworkers7 for large organic dye solutes. These first comparisons have been encouraging, in that reasonable agreement with experiment is obtained.19 By contrast, dielectric continuum-based theories3 are much less s u c c e s s f ~ l . ~ ~ ~ ~ ~ ~ In the present paper, we present the first molecular dynamics (MD) simulation test of the vdZ-H relation (eq 1.2). (For an earlier MD study of dielectric friction per se, see ref 21 .) The solutes selected for study are comparable in size to the polar solvent molecules, and this represents an exacting scrutiny (vide infra); by contrast, the experimental studies6*’have employed large probe molecules, for which the derivationlo of eq 1.2 is most plausible [cf. section 41. The outline of this paper is as follows. Theoretical background is provided in section 2, while section 3 details the model systems selected and the relevant M D procedures. Our results are presented and discussed in section 4, while section 5 concludes.
2. Theoretical Background A. Dielectric Friction. In the generalized Langevin equation (GLE), eq 1.3, for a rotating dipolar solute, the total time-dependent friction { ( t ) is related to the time correlation function of the so-called random torque exerted by the solvent on the solute:22
and the total friction constant { is its time integral
The dynamics of I?(t) are projection-operator-modified,I8 and hence are not standard Hamiltonian dynamics. As a result, { ( t ) and { are not accessible from direct simulation of the tcf in eq 2.1. However, {is related to the normalized angular velocity tcf ( 33(t ) )
C&) = (3.3)
(2.3)
(18) Mori, H. Prog. Theor. Phys. (Jpn.) 1965,33,423. Beme, 9. J.; Harp, G. D. Adv. Chem. Phys. 1970,22. Hynes, J. T.; Deutch, J. M. In Physical Chemistry: An Advanced Treatise;Eyring, H., Jost, W., Henderson, D., as.;
Academic Press: New York, 1975; Vol. 11.9, p 153. (19) Estimates of ground and excited state dipole moments are required in the implementation. See refs 6 and 10. (20) However, the situation for continuum theories improves when the finite character of the solute charge distributions is accounted for. See: Alavi, D. S.; Waldeck, D. H.J. Chem. Phys. 1991,94,6196; Hartman, R. S.; Alavi, D. S.; Waldeck, D. H. J. Phys. Chem. 1991, 95, 7872. (21) Cross, A. J.; Simon, J. D. J . Chem. Phys. 1987,86,7079. (After the present paper was submitted, we became aware of the work of Vijayakumar, P.;Tembe, 9. L. J . Phys. Chem. 1991, 95, 6430, in which the translational frictions on the members in an ion pair in a polar solvent are simulated, as are several solvation relaxation tcfs. The purposes of this reference are rather different than those of the present work.) (22) Note that a factor of ,3 = l/kTis missing in eq 3.12 of ref 10. The late Dr. T. Fonseca directed our attention to this misprint.
by23 i m d tC&) =
(2.4)
rl
and since the rotational diffusion constant Dratis also related to this time integral, we have23 Drot = kBT/I{
(2.5)
Thus, simulation of Ci(t)provides a route to { (and Drat).In addition, the time-dependent friction { ( t ) is also accessible from C;(t), via the Laplace transform relation +@) =
kJ+ ?@)I-lPt@)
(2.6)
which follows from eq 1.3 and in which @(t)
= -cz(t) =
(d. d ( t ))
(3.3)
and the Laplace transform is defined by, e.g.
Simulation ,Of @(t), Laplace transformation, and then solution of eq 2.6 for {(p) give {(r) by inverse Laplace t r a n s f o r m a t i ~ n . ~ ~ As will be further discussed in section 3, we will consider a neutral diatomic solute (which we call a neutral pair (NP)) as well as partially and fully ionic pair solutes. Our definition for the “dielectric” components of { and { ( t ) will be the difference from the neutral pair result, i.e.
Such definitions guarantee that the “dielectric” friction components vanish when there are no solute-solvent electrostatic effects, and do not require the assumptionz5 that nonelectrostatic and electrostatic tcf contributions to the friction are additive. B. Solvation Dynamics. The microscopic expression for the solvation time has already been given in eq 1.1. The explicit expression for the solvent variable 6AE is
6AE =
- (U)SP;
=
Vsolv.IP
- Vsolv.NP (2.10)
where ((...))sp denotes an equilibrium average in the presence of the particular solute pair (SP) and AE is, for a given solvent configuration, the difference in interaction energy Vsolv,sp due to all the solvent molecules between the fully ionic solute pair (IP) and a neutral solute pair (NP) [cf. section 31. It should be stressed that we use the same solvent variable, eq 2.10, for all three SP cases that we will consider. There are now many simulation studies of the 6AE(t) correlation function A(?), eq 1.1, for assorted solute and solvent systems.11-15,26-28 It is worth noting that there is also a GLE for the solvent variable GAE(t), in which a time-dependent friction on the solvent variable appear^."*^^-^^ This is not the same as the friction X t ) appearing in eq 1.3 and section 2.A; that friction refers instead to the solute nuclear coordinates. (23) See, e.g.: Boon, J. P.; Yip, S. Molecular Hydrodynamics; McGrawHill: New York, 1980. (24) See,e.g.: Karim, 0. A.; Haymet, A. D. J. J . Chem. Phys. 1990,93, 5961, for a related example. For a method to accomplish the inversion, see: Zichi, D. A.; Kim, H. J.; Carter, E. A.; Hynes, J. T. To be submitted. (25) See, e.g., ref 6. (26) Carter, E. A.; Hynes, J. T. J . Phys. Chem. 1989, 93, 2184. (27) Zichi, D. A.; Ciccotti, G.; Hynes, J. T.; Ferrario, M. J. Phvs. Chem. 1989, 93, 6261. (28) (a) Smith, B. 9.; Hynes, J. T. To be submitted for publication. (b)
Smith, B. B.;Hynes, J. T. Proceedings of the Trieste Workshop on Electrochemistry; in press. (c) Smith, B. B.; Kim, H. J.; Borgis, D.; Hynes, J. T. Mataga Conference on Elecrron Transfer; in press. (29) Hynes, J. T. J . Phys. Chem. 1986, 90, 3701. (30) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. J . Chem. Phys. 1990, 93, 7137.
Bruehl and Hynes
4070 The Journal of Physical Chemistry, Vol. 96, No. 10, I992
NP HIP IP
0.31 f 0.03 0.44 f 0.05 0.58 f 0.1
14.8 f 0.8
33.5 f 6.7 105.0 f 19.0
0.078 f 0.004 0.035 f 0.007 0.011 & 0.002
95.8 f 0.5 178.0 i 3.6 361.0 f 6.5
18.7 & 3.8 90.2 f 16.0
82.2 f 3.0 265.2 f 17.2
C. vdz-H Connection. With the above definitions, the vdZ-H connection between the solvation time and the dielectric friction assumes the form 7M = sg/sy(t=o) (2.1 1) the components of which can be calculated by simulation, as described next. 3. Solute and Solvent Models and Simulation Method Three equilibrium, constant-temperature3I molecular dynamics (MD) simulations are performed for a solute pair immersed in a dipolar solvent. The solvent28*30~32~33 is composed of 342 rigid dipolar diatoms constructed from Lennard-Jones (LJ) spheres of mass 40 amu separated by 2.0 A and containing point charges of f e / 4 such that the dipole moment is 2.4 D. At 250 K and number density 0.012 A-3, this solvent is meant to roughly approximate methyl chloride, a representative polar aprotic solvent. The three MD simulations differ in the dipole moment of the SP. Each SP is constructed from LJ spheres of mass 40 amu rigidly separated by 3.0 A. We simulate an ion pair (IP) with point charges of f e making p = 14.4 D, a one-half ion pair (HIP) with point charges of f e / 2 making p = 7.2 D, and a neutral pair (NP) with no charges. The total potential energy is a sum of LJ pair interactions and Coulomb potentials V(rij)= VLJ(rij) + zizjrij-I between particles of charge zi and z, separated by rip The LJ spheres making up the components of both the solute and solvent diatomic molecules are described by LJ parameters e/ke = 200 K and a = 3.5 A. The simulation is carried out in a cubic box of length 30.52 8, whose boundary conditions are periodic. We integrate the equations of motion according to the Verlet algorithm34with a time step of 0.01 ps. The SHAKEalgorithm35is used to maintain the bond length constraints. The long-range forces are treated using the Ewald summation method.36 At the start of each simulation, the equations of motion are integrated for 20 ps (2000 time steps) in order to verify that equilibrium has been attained. After this initial period, the simulation is continued for a total of 500 ps, during which the quantities needed for the time correlation functions are collected. These include the angular velocity, angular acceleration, and dipole direction of the SP and the solvent response to the rotating SP, namely, the value of AE. When the simulation is completed, the time correlation function of the fluctuations of each of the quantities is constructed according to standard methods. It should be noted that the angular yelocities and accelerations arecalculated from 3 = ti X ti and 3 = ti X ti. The quantities ti, ti, and ti are the dipole direction, velocity, and acceleration vectors which are resolved in a body-fixed, nonrotating coordinate system. In order to keep constant the moment of inertia of the SP, the angular velocities and accelerations must be resolved in a body-fixed, rotating frame whose origin corresponds with the SP center of mass and whose axes fall along the principle axes of r~tation.~'Once this is accomplished, the moment of inertia tensor is diagonal and independent of time, with two of the components equal and the third component, which is coincident with the molecular axis, zero. (31) Nose, S. J . Chem. Phys. 1984, 81, 511.
(32) Borgis, D.; Hynes, J. T. J. Chem. Phys. 1991, 94, 3619. (33) Ciccotti, G.; Ferrario, M.; Hynes, J. T.; Kapral, R. Chem. Phys. 1989, 129, 241.
(34) Verlet, L. Phys. Rev. 1967, 159, 98. (35) Ciccotti, G.; Ryckaert, J. P. Comput. Phys. Rep. 1986, 4 , 345. (36) Hansen, J. P. In Molecular Dynamics Simulations of Statistical Mechanical Systems; Ciccotti, G., Hoover, W. G., U s . ; North-Holland: New York, 1986. (37) Goldstein, H. Classical Mechanics, 2nd ed.; Addison-Wesley: Reading, MA, 1980.
050
0.25
iI
' &' I6
0
I
1
2
3
4
5
1. Solvation tcf A(?) vs time in ps for IP (-), HIP (---),and NP The initial Gaussian eq 4.1 contributes 2996, 452,and 66% of the total time areas, respectively. The long-time tails of the tcfs are well tit by the following exponential functions: 0.24 exp(4.521) (IP); 0.30 exp(-O.95?) (HIP); 6.0 exp(4.9) (NP). For clarity, representative estimated error bars3*are floated above the curves. F i i (.a*).
The solvation tcPs in eq 1.1 have been previously calculated for the H I P in ref 27 and for the N P and IP in ref 1 la. They have been recalculated in the present work during the same simulations that generate the frictional quantities.
4. Results and Discussion A. Solvation Dynamics. The results for the normalized 6AE tcf, A(?), are displayed in Figure 1. There are several features to note. First, the significant initial drop in A(t) is very similar for the NP, HIP, and IP solutes." This initial decay is well described by a Gaussian time dependence]' &(t) = exp(-wsp2t2/2)
(4.1)
where wsp is the solvent frequency (4.2) As described in ref 11, this Gaussian behavior is related to the inertial free streaming of the solvent molecules, whose ensuing coordinate changes alter AE. These solvent translational and orientational changes are independent of the interaction potentials. Such Gaussian behavior has now been found in a number of s t ~ d i e s . ~ For ~ - ' ~the current solute pairs, the solvent frequency has the values, in ps-', of 6.1 (NP), 6.5 (HIP), and 7.7 (IP). The relative insensitivity of wsp to the solute identity is discussed in refs 11 and 26 in terms of its relationship via wsp = (ksp/msp)'/2 to the solvent force constant and mass. The longer time behavior of A(r) depends on the intermolecular interactions" and differs noticeably for the various SP cases, with the slowest relaxation for the full IP solute. It is this nonlinear feature that lengthens the solvent relaxation time rm [cf. Table I] as the dipole moment of the SP increases." B. Rotatio~lDynamics. Figure 2 displays the normalized angular velocity tcf C,(r), eq 2.3, for the three SP cases. There is a weak caging behavior for the N P the angular velocity at time t is negatively correlated to its initial value. This correlated reversal of the angular velocity direction is related to the nonspherical (38) The estimated error bars for the tcfs are calculated according to: Zwanzig. R.; Ailawadi, N. K. Phys. Reu. 1969, 182, 280.
The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 4071
Dielectric Friction and Solvation Dynamics 1.o
0.5
C,( t ) \
0
\
-0.5 I
0
0.5
1.5
1
2
t (PS) Normalized angular velocity tcf C&) vs time in ps for IP (-), HIP (---),and N P For clarity, representative estimated error bars are floated above the curves. F v 2.
( e ) .
2
0
8
6
4
10
Figure 4. Orientational tcf C2(r)vs time in ps for the IP (-), HIP (---), and NP (-). The orientational t c f s are fit by the following exponential functions: exp(4.076t) (IP); exp(4.27t) (HIP); exp(4.59r) (NP). Representative estimated error bars are included on the IP curve.
I
i
.,..,..---& - \
0.0'
0
"
2
"
4
"
6
"
a
" 10
t (PSI Figure 3. Orientational tcf C,(r) vs time in ps for IP (-), HIP (- - -), and N P (-). The orientational t c f s have been fit to the following exponential functions: exp(4.025t) (IP); exp(4.077t) (HIP); exp(4.19t) (NP). Representative estimated error bars are included on the IP curve.
geometry of the NP and is not dissimilar to that observed in simulations of nonpolar linear molecules.39 However, as the dipolar character of the SP increases, the electrostatic interactions with the model dipolar aprotic solvent molecules cause a much stronger rotational caging, reaching a pronounced oscillatory pattern for the fully charged IP case. Such strong caging behavior has been observed in simulations for various polar and/or charged solutes in quite different solvent systems, e.g., hydrogen-bonding solventsa and molten salts.41 These features result in a marked increase in the rotational friction constant t [cf. eq 2.41-and a corresponding decrease in D,,, via eq 2.5-as Table I indicates. The defined dielectric component SY = f - tNp of the friction is already comparable to the NP friction tNp for the HIP case, and exceeds it by a factor of about 6 for the IP case. For the IP case, this dominance of (39) Ruhman, S.;Nelson, K. A. J. Chem. Phys. 1991,94,859. Quentrec, B.; Brot, C. Phys. Rev. A 1975, 12, 212. Berne, B.J. In Physical Chemistry: An Advanced Treatise; Eyring, H., Jost, W., Henderson, D., Eds.; Academic Press: New York, 1975; Vol 8.B, p 539. (40) (a) Karim, 0. A.; Haymet, A. D. J. J . Chem. Phys. 1990,93,5961. (b) Matsumoto, M.;Gubbins, K. E. J. Chem. Phys. 1990, 93, 1981. (c) Frattini. R.; Ricci, M. A.; R u m , G.; Sampoli, M. J. Chem. Phys. 1990, 92, 2540. (d) Evans, M. W. J. Chem. Phys. 1987,86,4096. (e) Wojcik, M.; Clementi, E. J. Chem. Phys. 1986,85, 3544. (f') Spohr, E.; Heinzinger, K. J. Chem. Phys. 1986,84, 2304. (41) (a) Okaulki, S.;Ohtori, N.; Okada,I. J. Chem. Phys. 1990,93,5954. (b) Kato, T.; Machida, K.; Oobatake, M.; Hayashi, S.J. Chem. Phys. 1990, 92, 5506; J. Chem. Phys. 1988,89,1411.
0
-
0.5
1
t
1.5
2
(PS)
Figure 5. Time-dependent friction { ( r ) in PS-~vs time in p for IP (-),
-
HIP (- -), and N P (-). {(r) is calculated as described in section 2.A. Though not emphasized in the text, the initial decays here are well described by Gaussian functions, which contribute 5 5 1 , 59%, and 79% to the total friction constant. For clarity, representative estimated error bars are floated above the curves.
the total friction by its dielectric component is a contribution considerably greater than that estimated for large organic dye solutes.6 We have also determined D,,,and t by an alternate route as a check. In particular, we have calculated the MD orientational correlation functions C,(O = (P.P(o) = ( P J ~ w~ i S) ; C2(t)= W ( P . P ( ~ ) ~11) = (p2[cosWl)
(4.3)
in which P.P(r) = cos e ( t ) is the cosine of the angle between the solute dipolar axis at time t and its initial value. If the solute reorientational motion were d i f f ~ s i o n a l , "then ~ ~ ~ C,(t) ~ and C2(t) would decay e~ponentially~~
C,(t) = exp[-l(l + l)Drott] (4.4) Figures 3 and 4 show that, after a brief transient period, exponential behavior is in fact well satisfied. This is itself noteworthy, since this behavior is not always observed for solute rotation.aa The rotational diffusion constants defined by the application of eq 4.4 are found to be the following: NP, 0.104 ps-I (Cl), 0.088 (42) Fixman, M.J. Chem. Phys. 1968, 48, 223. (43) See,e&: Fleming, G. R. Chemical Applications of Ultrafmt Spectroscopy; Oxford: New York, 1986; FrBlich, H. Theory of Dielectrics; Clarendon: Oxford, 1958.
4072 The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 300
Bruehl and Hynes
1 -
1.0
R
HIP
@ j j
200
I\
I
tcf 0501
!\
0
1
1
\
( C ---__
\
0
05
1
t (PSI
2
0.0
1
0
15
0.5
\
1.5
1
2
Figure 6. 6 { ( f ) in P S - ~vs time in ps for IP (-) and HIP (- -). &{(to is the dielectric component of the friction as defined by {(f) - tNp(f).For
t (PS) Figure 7. Solvation tcf A(f) (-) and normalized time dependent dielectric friction 6f(f)/&{(f=O) (-) vs time in ps for HIP. For clarity,
clarity, representative estimated error bars are floated above the curves.
representative estimated error bars are floated above the curves.
-
TABLE II: Solvent Relaxation-Dielectric Friction Comparison HIP
IP
T A F . DS 0.44 f 0.05 0.58 f 0.10
6r/&t(f=o), DS 0.23 f 0.04 0.34 f 0.05
ps-I (C2);HIP, 0.036 ps-' ( C ] ) 0.035 , ps-I (C2);NP, 0.013 ps-I ( C , ) ,0.012 ps-I (C2). These compare well with the D,,, values obtained from the CJt) time area [cf. Table I]. Finally, the time-dependent-friction values, eq 2.1, calculated as described in section 2.A, are displayed in Figure 5, and the corresponding dielectric components a{ = { ( t ) - lNP(t) for the dipolar solutes are shown in Figure 6. Comparison of Figures 5 and 6 with, e.g., Figure 3 for C,(t)indicates there is no important reorientational motion on the time scale of the lifetime of the various frictions. It is noteworthy that the initial value {(t=O) of the friction is responsible for the major portion of the variation of {in the NP-IP sequence. If we define [cf. eq 2-21 { = {(f=O)~r
(4.5)
then Figure 5 and Table I show that while the time area 7f increases by slightly less than a factor of 2, {(t=O) increases by a factor of almost 4.44945A similar feature holds for the dielectric component: on going from the H I P to the IP, the initial value s{(t=O) slightly more than triples, while the time in
a{
= 6{(f=0)76r
(4.6)
is increased by a factor of 1.5. C. vdZ-H Relation Test. With all the required quantities in hand, we can now test the molecular vdZ-H connection (eq 2.1 l), and Table I1 details the results. On the positive side, the vdZ-H relation is satisfied to within a factor of 2 and also properly accounts for the approximately 50% lengthening of both 7M and resulting from the doubling of the dipole moment of the HIP to obtain the IP. The latter feature would be completely missing (44) The analogue of this initial friction in chemical reaction problems is central in determining the deviation of the reaction rate constant from its transition state theory approximation for sharp barrier reactions [refs 16a.b and 171. (45) The product of this initial friction and the moment of inertia I is suggested in ref 10 to be equal to an experimental relaxed fluorescence red shift. The latter can be extracted from the simulation results of ref 1l a [from the Stokes s h i f t s relation: s = (1/2)(k1p + k ~ p ) ( ( A i ? ) ~ p- (hE)Np)*] as 19928 cm-I for the NP-IP fluorescence system, while 16{lp(t=O) from the present work is 39912 cm-I. For the HIP case, if we estimate the force constant kHIpby the average of klP (8.0 X lo-' cm-') and kNP (4.9 X 10'' cm-I), then with ( A E E ) H= ~ ~9531 cm-' and ( A E ) N P= 71 cm-I from the present work, the shift is estimated as S = 5101 cm-I, while 16fHlp(t=O) = 12 370 cm-'. These estimates indicate that the Stokes shift will underestimate the initial dielectric friction value by about a factor of 2, at least for comparably sized solute and solvent molecules.
It
I
0.0 0
0.5
1.5
1
2
t (PSI Figure 8. Solvation tcf A(f) (-) and normalized time dependent dielectric friction 6 { ( f ) / 6 { ( f = O ) (-) vs time in ps for IP. For clarity, representative estimated error bars are floated above the curves.
in a dielectric continuum a p p r o a ~ h , ~where - ' ~ all relaxation times would be independent of the solute dipole moment. On the negative side, the factor of 2 discrepancy in Table I1 is nonnegligible. We can pursue this somewhat further by consideration of the behavior of the time integrands which enter into eq 2.11; Le., we consider the validity of the relation given by VdZ-H" (4.7)
Figures 7 and 8 show that while there is reasonable coincidence of the two functions for short times, there is significant deviation at longer times f 2 0.4 ps. In particular, the longer time tail of A(f) is missed by the normalized timedependent dielectric friction, and this feature is largely responsible for the underestimate of r M by 7af [cf. Table 111. The 'missing" tail in 8{(t)-which is not the result of a cancellation between s,,(~) for the half- and the full-ion pairs and the neutral pair lNP(t),cf. Figures 5 and 6-is quite striking from a dielectric continuum point of view. If one imagines that a continuum description is most likely valid at longer times, then one would expect an equal amplitude tail for each of A(?) and 6f(t):lO
A(t)
-
8f(f)/6f(t=O)
-
exp(-t/rL);
7~ =
2t,
+1
2to + 1 lD
(4.8)
Here 7Lis the solvent longitudinal relaxation time, rDis the Debye
The Journal of Physical Chemistry, Vol. 96, No. 10, 1992 4073
Dielectric Friction and Solvation Dynamics time, and e, and eo are the solvent high frequency and static dielectric constants, respectively. This common asymptotic behavior is clearly not applicable to the current comparably sized solute and solvent molecule system.& A further point of interest is that the relatively slow orientational relaxation of the HIP and IP solute pairs appears not to be involved in the long-time solvent relaxation exhibited by A(t). For example, the tail of A(t) for the IP decays approximately 7 times faster than Cz(t)[cf. Figures 4 and 81. Some insight on the more rapid decay of 6{(t)/6{(t=O) compared to that of A(t) and the difficulties with the vdZ-H relations eqs 2.1 1 and 4.7 may be provided by considering the limit of a point dipole solute. In this case, the random torque in the frictional tcfs is of the form fl/*Bd(l)= i i ( r ) x E ( t )
it
x i(t)
(4.9)
where E ( t ) is a random electric field experienced by the solute. In eq 4.9 the slower solute reorientational dynamics is neglected and projection-operator-modified (POM) dynamicsIs is unierstood. The point dipole approximation gives M(r)= -ii.E(t), again ignoring the rotation and where POM dynamics are not involved. In this case then-aside from POM dynaFics issues-the dynamics of both tcfs are governed by those of E(t), as in the vdZ-H treatment.I0 For a finite dipole, however, this is not so. While AE involves the electric potential at the dipolar sites [cf. eq 2.101, the random torque will involve a spatial gradient of the potential, Le., a random electric field, at each site. It is then not unreasonable to expect that the random torque correlation could decay more rapidly than that of 6AE. The POM dynamics issue is more difficult to address. Their explicit presence in 6 K t ) but not in A(r) is ignoredl0 in eqs 2.1 1 and 4.7. If, for example, the random torque decayed rapidly compared to the angular velocity correlations, the consequences of this difference would be expected to be minima1.18*47For the present comparably sized solutesolvent systems, with pronounced caging, comparison of Figures 2 and 5 shows that this assumption is not satisfied. For solutes larger and more massive than the solvent molecules, this difficultyshould subside (although the finite charge distribution aspects discussed above could remain); ultimately the friction would reduce to that for a fmed solute. Simulations for larger solutes would be of interest in connection with all these issues. In summary, the vdZ-H relations eqs 2.1 1 and 4.7 implicitly assume a point dipolar charge distribution for the solute and a generalized Brownian motion neglect of POM dynamical motion of the solute. In view of the above discussion, it is useful for perspective to note that the longer time discrepancies between A(t) and S{(t)/6Rt=O) can often be largely irrelevant in a chemical reactionlo context. To illustrate this,we display in Figure 9 the Grote-Hynes transmission ~ o e f f i c i e n t ~ ~ J ~ ~ ~ ~ * ~ * K
= [K
+ Ub-IJmdf
exp(-UbKf) { ( f ) ] '
(4.10)
for a model parabolic barrier reaction (with barrier frequency 4). The reaction system is assumed for the present illustrative purposes to be subject to the same friction as experienced by a rotating ion pair. The K values are separately calculated using both sides of eq 4.7, i.e., for i = HIP, IP (4.1 1) li(f)= {NP(~) + S{i(t=o)[6{i(t)/6S;(t=o)l {i,A(t)
= {NP(r)
+ a{i(t=o)
Adt)
(4.12)
The former friction is the correct one within the model considered follows from the approximate here, while the second friction c,;;a(t) equation, eq 4.7. Figure 9 shows that, over most of the barrier frequency range, there is quite good agreement using the approximate {i,A. It is (46) For a related failure of simple asymptotic time behavior, see ref 1 1 . (47) See, e.g., ref 17a and references cited therein. The validity of the vdZ-H relations in the point dipole, fixed solute limit has now been considered by B. Bagchi (1992 preprint). (48) Grote, R. F.; Hynes, J. T. J . Chem. Phys. 1980, 73, 2715.
a) HIP
I
n
I
0
20
40
60
80
[ a+,/(kBT/l) "'1 Figure 9. Transmission coefficients K versus reduced barrier frequency ~ ~ / ( k ~ T /for r ) (a) l/~ HIP and (b) IP cases. (-) KGH from the GroteHynes (GH) equation (eq 4.10), with the actual model friction equation (eq 4.1 1); (---) KGH with the approximation friction equation (eq 4.12); (-) Kramers transmission coefficient equation (eq 4.13) with the friction constant asscciated with the actual model friction equation (eq 4.1 1); (-) Q~ calculated with the approximate friction equation (eq 4.12), except that the prefactor initial value b(,(t=O) is approximated by a vdZ-H relationlodescribed in ref 45. At the lowest barrier frequencies [wb/ (knT/r)*/*d 3.8 for the HIP and ;5 2.0 for the IP], the GH and Kramers results coincide to within 15%. In that frequency range, K becomes m a t sensitive to errors in the friction constant.
only at the lowest barrier frequencies that the transmission coefficient becomes sensitive to the longer time tail present in A(t) but not in cSy(t)/SR.(t=O).The extent of the sensitivity of K to the short-time behavior of the time-dependent friction is further highlighted in Figure 9 via the Kramers theory transmission coefficient49
which is only sensitive to the full time area of the friction, i.e., the friction constant. 5. Concluding Remarks We have shown, by molecular dynamics simulations, that the van der Zwan-Hynes relationship connecting the solvation time and the dielectric friction holds to within a factor of 2 for halfand full-ion pair solutes in a model dipolar aprotic solvent. In addition, the trend that the solvation times and dielectric friction constants increase in concert with more dipolar solutes is correctly reproduced. The deficiencies of the relation for this stringent case of comparably sized solute and solvent molecules are traced to a longer time tail in the solvation dynamics that is absent in the timedependent dielectric friction. This tail, which is presumably collective in character, awaits a detailed molecular explanation. (49) Kramers, H. A. Physica 1940, 7, 284.
J. Phys. Chem. 1992, 96,4074-4076
4074
Its marked dependence on the solute pair identity indicates to us that a successful theory must explicitly incorporate microscopic details of the solutesolvent interaction. Similar remarks apply to the time-dependent friction itself. Acknowledgment. This work was supported in part by N S F
Grant CHE88-07852 and by a grant of Cray YMP time from the Pittsburgh Supercomputing Center. M.B. thanks the N S F Postdoctoral Fellowships in Chemistry Program and notes that this material is based upon work supported by the National Science Foundation under Grant CHE-9002183 awarded in 1990. We thank Professor John Simon (UCSD) for useful discussions.
From the Density of States to the Velocity Correlation Function in Liquids G. Seeley, T. Keyes,* and B. Madan Department of Chemistry, Boston University, Boston, Massachusetts 0221 5 (Received: November 7, 1991)
Harmonic normal mode analysis (NMA) is used to estimate the Fourier transform, denotedfiw), of the velocity correlation function in liquids. Increasingly good agreement betweenfiw) and the NMA density of states is obtained with increasing pressure, in agreement with the contention that NMA is a powerful technique for liquids with low fluidity. The work is in the spirit of Fixman’s effective harmonic approximation, in which NMA is applied to systems in which it is not literally valid.
Introduction The harmonic approximation, with its associated normal mode analysis (NMA), provides a complete description of both the static and dynamic properties of low-temperature crystals. We will discuss atomic systems with N atoms. The total potential energy, U,is expanded U(r) = U(ro)+ F(ro).x
+ f/zC(ro):xx
(1)
where r denotes the 3 N atomic coordinates, ro the lattice coordinates, and x the displacements of the atoms from the sites, r = ro x. The lattice configuration is a minimum of U, so F(ro) = 0 and C is the force constant matrix. As T 0 the atoms do indeed stay close to the lattice sites, x is the small, and the truncation at second order is literally correct. The eigenvalues of G are the squares of the natural frequencies of the system, and binning the frequencies yields the density of vibrational states, p(w). The eigenfunctions are the normal modes of vibration, e.g., the phonons. The simplicity of this description suggests that one try to extend it to situations in which it does not obviously hold. Fixman’ has introduced, with considerable success, an ‘effective harmonic approximation”, in which a best harmonic representation of a finitetemperature anharmonic crystal, one in which the truncation at second order is not strictly valid, is constructed. This work was extended by Wolynes et al.2 In a similar spirit, we have p r o p o ~ e d that ~ . ~ harmonic N M A be applied to liquids. One might think this impossible due to the absence of lattice sites, but such is not the case. In this article, we will show that N M A is a promising route to the velocity correlation function in liquids.
-
+
Harmonic Normal Mode Analysis in Liquids
We begin by introducing the harmonic density of states for a liquid. In the absence of lattice sites we evaluate the force constant (second derivative of U) matrix C ( r ) for an arbitrary liquid configuration r, obtain p(w;r), and then perform a standard configuration average, denoted by angle brackets, to obtain ( ~ ( w ) ) . So far we have calculated ( p ( w ) ) via computer simulation and numerical eigenanalysis, although Stratt5 has given a method (1) Fixman, M. J . Chem. Phys. 1969, 51, 3270. (2) Stoessel, J. P.; Wolynes, P. G. J . Chem. Phys. 1984, 80, 4502. (3) Seeley, G.; Keyes, T. J . Chem. Phys. 1989, 91, 5581. Madan, B.; Keyes, T.; Seeley, G. J . Chem. Phys. 1990, 92, 7565; 1991, 94,6762. (4) Seeley, G.; Keyes, T.; Madan, B. J . Chem. Phys. 1991, 95, 3847. (5) Xu, B.; Stratt, R. J . Chem. Phys. 1990, 92, 1923.
0022-36S4/92/2096-4074$03.00/0
which could lead to accurate nonsimulational approximations. Since r is unlikely to be a minimum of U,C ( r ) has negative eigenvalues and imaginary frequencies. If the harmonic expansion in liquids were to be taken literally, which it is not, real frequencies would correspond to stable modes and imaginary frequencies to unstable modes. We plot the contribution of the imaginary frequencies to ( p ( w ) ) along the negative real axis. The density of states vanishes a t zero frequency, and thus ( p ( w ) ) has distinct stable and unstable mode lobes, ( p , ( o ) ) and (pu(w)). For ( p ( w ) ) normalized to unity, the area of (p,(w)) is the fraction of unstable modes, denotedf,. We have shown3” thatf, is a sensitive indicator of fluidity, and we have given theories of the self-diffusion coefficient, D, in the supercooled regime, in whichf, is the crucial T-dependent3 or P-dependent4 parameter. In this paper we wish to use NMA to continue the study of pressure dependence of liquid-state dynamics at constant T begun in ref 4; our earlier work3 studied the T-dependence at constant density. In Figure 1 we show ( p ( w ) ) at 300 K and three pressures in Lennard-Jones argon. The decrease in fu as P increases and the system becomes less fluid is evident; note also the shift of ( p , ( w ) ) to higher frequencies. What can we do with ( p ( w ) ) ? A route to physical dynamical quantities is suggested in a paper by Zwanzig6 and in the work of Stillinger.’ These latter authors considered the motion of the point representing all 3 N coordinates of an N-particle atomic fluid (phase point) in the 3N-dimensional configuration space, At low T the point spends long,periods in the vicinity of a local minimum of U,quickly crosses a barrier to a neighbor minimum, and so on. Stillinger et al. formulated and verified this picture in simulations by continuous quenching, in which one always knows the energy/particle, c#J(~), of the local minimum to which the phase point will drain. They verified that the barrier hopping rate, denoted wh, does indeed decrease sharply as Tis decreased. The P-dependence of the hopping rate is yet to be reported. The most important dynamical quantity in atomic systems is the velocity correlation function, written kTC(t), whose time integral equals D. In a true harmonic crystal, the Fourier transform of C ( t ) , denoted f i w ) , equals p ( w ) . In a harmonic amorphous system with an ensemble of local minima, flu) = ( ~ ( w ) ) where , the q denotes an average over the minima. Zwanzig suggested5that barrier hopping produces decorrelation of the contributions of the harmonic modes to C(t). The resulting formula is just the harmonic C ( t ) multiplied by a “lifetime (6) Zwanzig, R. J . Chem. Phys. 1983, 79, 4507. (7) Stillinger. F. H. J . Chem. Phys. 1988,88,7818 and references therein.
0 1992 American Chemical Society
