J. Phys. Chem. 1993, 97, 13-17
13
A Simple Interpretation of Polar Solvation Dynamics M. Maroncelli,' Vijaya P. Kumar, and Arno Papazyan Department of Chemistry, 152 Davey Laboratory, The Pennsylvania State University, University Park, Pennsylvania 16802 Received: October 1 , 1992; In Final Form: November 17, 1992
In this letter we describe a n approximate relationship between the solvation response to changes in solute charge, C v ( t ) ,and the reorientational dynamics of molecules in the pure solvent, C l ( t ) . Computer simulations show that the function {Cl(t)),s,where asis the solvent dipole density, provides a good approximation to C v ( t ) in a wide variety of systems. This result leads to a simple and intuitive interpretation of the solvation response in terms of molecular motions.
in terms of the dipole
Introduction The phenomenon of polar solvation dynamics, the response of a polar solvent to changes in the charge distribution of a solute, has been the subject of many experimental and theoretical studies over the last decade.'-4 Two of the most characteristic features of the solvation response are its surprising rapidity and the fact that the rate of this response is intimately linked to the polarity of the solvent. Even though the time dependence of solvation mainly reflects reorientational motions of solvent molecules, the time scale for solvation energy relaxation is much faster than times normally associated with single-particle reorientations. The greater the polarity of the solvent, the more "collective" is the solvation response and the faster it is relative to single-particle dynamics. In this letter we suggest an approximate formula that captures this polarity-dependent connection between the solvation response and the rotational dynamics of solvent molecules in a simple and intuitively appealing way.
-
The dynamic we consider is the response to a step-function change in the charge (Qi Qf)of a stationary monatomic solute. The time dependence of this response is described by the normalized function
where ( V(t)) is the average reaction potential (electrical potential as the solute charge site) at time t . Rather than model this nonequilibrium property directly, it is more convenient to shift attention to the equilibrium time correlation function (tcf'):
which is equivalent to S,-r(t) in the limit of a linear solvent response.' In this expression, bVis the fluctuation in the reaction potential, V - ( V),, and (X),denotes the average value of the property X in a system in equilibrium with a solute of charge a = i or f. Computer simulations have shown that although the linear response assumption is not exact, Cv(r) provides a very We good approximation to Si-.r(t) in many cases of therefore restrict our attention to the linear response limit and seek to understand the solvation response by modeling this solvation correlation function. In particular, we relate this collective tcf to the reorientational dynamics of individual solvent molecules in the absence of the solute. The latter are specified
(p)
autocorrelation function
(3) The Solvation Frequency We begin by examining Cv(t)at early times where thedynamics results from independent "free-streaming" of solvent molec~les.~3~ What is meant by this terminology is the following. At any given time solvent molecules possess velocities representative of an equilibrium Maxwell-Boltzman distribution. In a classical system this distribution is completely independent of intermolecular interactions. At sufficiently short times, the decay of thesolvation tcf reflects the average effect that free propagation of these random velocities has on the reaction potential. When this free-streaming prevails Cv(t) is found to be well approximated by a Gaussian function, Cv(t)
= exp(-'/2ws 2t 2)
(4) where ws is the "solvation frequency", or the t = 0 curvature of C"(t) (SV22,
with = dV/dt.S,7 One of the surprising results of computer simulations on small-molecule solvents has been the observation that this free-streaming phase of the motion often dominates the decay of Cv(t). In most simulations, solvation tcfs have been observed in which the initial Gaussian phase accounts for well over 50% of the total r e l a x a t i ~ n . ~Given - ~ ~ this ~ ~ ~fact, ~ a model capable of predicting the solvation frequency would be a useful first step toward understanding the solvation response. A general expression for the solvation frequency can be simply derived if one is willing to accept two important assumptions. First, weassume that the translational contribution to thesolvation dynamics can be ignored. While translational motions of solvent molecules near to the solute are expected to play some role in the solvation response,'-3 simulations s ~ g g e s t ~that * ~inJ ~many cases such motions are not of primary importance. At least in the systems that have been studied in detail, the motions occurring during the free-streaming period have been found to be primarily reorientational motions identical to those that would be observed in a gas of noninteracting solvent molecules.s*'O Thus, for the sake of simplicity we consider only the reorientational dynamics. Second, we also ignore the details of the charge distribution of solvent molecules and assume that the effect of solvent motions on the reaction potential can be modeled in terms of the rotation of point dipoles, U n d e r these two assumptions t h e
0022-365415812097-0013S04.00/0 Q 1993 American Chemical Society
14
Letters
The Journal of Physical Chemistry, Vol. 97, No. 1, 1993
TABLE I: Solvation Frequencies u,and Free-Rotor Frequencies 01 Observed in Computer Simulations
reaction potential can be written as Pj'rj "ri V=C;=pC, i ri i ri
(6)
CH,CN
where ri is the vector between the solute and solvent molecule i , and uri is the radial orientation of solvent molecule i, defined by uri = (p,-ri)/(plril). For ease of notation we consider fluctuations about a neutral solute for which (V) = 0. In this case (6V) is just
(7)
In the absence of translations the (6p) term in eq 5 takes on an identical form except that ti,; = du,i/dt replaces u,. However, the independence of coordinates and velocities means that, upon averaging, all i # j terms vanish from the double summation. (6p) can be thus related to the angular velocity of a single solvent molecule, w l , via the relations
In these expressions 81 denotes the time derivative of the unit vector along the dipole direction of molecule 1, and 01 is defined by u~ = wl X UI. The term in brackets in the second equation is equal to (6V),,,,the value of (SV) in a reference solvent having the same translational structure as the real solvent but in which the orientations of molecules are uncorrelated. Combining eqs 5 and 8, the solvation frequency can therefore be written as2=
WIC
solvent
(9)
with
Equation IO is one of the key results of the present work. It states that the initial decay of the solvation energy is related to the free rotational motion of solvent molecules ((01~)) by a translation factor cy,. The frequency ( w 1 2 ) is the average squared rotational frequency in a direction perpendicular to the dipole direction and it is determined entirely by the inertial properties of the solvent molecules. For example, in h e a r molecules ( w 1 2 ) = 2keT/I, where I is the moment of inertia and ksT is Boltzmann constant times the temperature. The factor cy, that translates between single-particle motion and the solvation response is the ratio of electrostatic fluctuations in the real system to those that would occur if solvent orientations were uncorrelated. That is, a, measures the extent to which orientational correlations serve to quench fluctuations in the solvation energy. An explicit expression for a,can be obtained by resorting to continuum solvent models for calculating (6V) and (bV),,.We assume that outside of some solute radius 'a" the solvent is a structureless continuum fluid with dielectric constant €0. The usual Born treatment of solvation then gives (6V) asll
Furthermore, for a continuum solvent the structural average appearing in (SV),, is simply
where p is the solvent number density. The factor cys can thus
solute"
SO
ref*
(ps-1) 3.14
1
S+
H20 'CHICI"
so, LO, L+ NP IP
Stockmayer +le Stockmayer +le DD O.le SO DD O.2e DD 0.3e DD 0.4e DD O.5e
14.0 16.4
17.1
S-
2 3
-I00
-23
4b
3.45 3.45 4.45
5
1.76
4a
6.2
2.28
=
a,
WI
(ps-')
7.4 6.9
8.3 -14
2.4 3.9 5.4 7.4 9.5
4wpp2/3kgT 19.9 19.7 27.3
w,~/wI*
29.7 -20 7.4
10.5 4.0 5.8 -9.4 1.9 4.9 9.3 18 29
18.6 8.39 4.82 4.82 8.99 (1.26) 5.02 11.3
20.1 31.4 Except for theCH3Clsimulations (see text) all solutes are monatomic. So,S+,and S-denote neutral and singly charged solutes. In water 'small" (S) uncharged solutes and all large (L) solutes showed approximately the same Cv(t) dynamices6 The solutes designated +le indicate nonequilibrium simulations in which the response S,l,(t) to charging a neutral atom was monitored. (See the original references for details.) 1 = ref 5, 2 = ref 6, 3 = ref 7, 4a = ref 9, 4b = ref 10. The entires labeled 'DD" are from previously unpublished simulationsof the 'dipolar dumbbell" solvents described in the text. In all cases except water w1 is the root-mean-squarefree rotational frequency w1 = (2kgT/fl1/*.For water 01 is the average rotational frequency about the two axes not coincident with the dipole direction. be written
Equation 13 is the second main result of this work. It expresses the fact that although the strength of solute-solvent interactions or (SV) would scale with the dipole density (term in brackets) if solvent molecules did not interact, the solvation energy is in fact much more weakly dependent on solvent polarity (i.e., grows only as 1 - e0-I) due to 'screening" by solvent-solvent interactions. For strongly polar solvents (€0 > 10) the second term in eq 13 is close to unity and the translation factor a,is essentially the dipole density ( 4 7 r p ~ * / 3 k ~ Titself. ) It should be noted that eq 13 should be more accurate than implied by its derivation from continuum approaches. Since only the ratio (SV),,/(6V) is of interest here, the well-known inadequacies of a continuum approach are less serious than in calculating the solvation energy itself, as suggested by the disappearance of the solute radius in the final result. The accuracy of eqs 9 and 13 can be tested against the results of molecular dynamics simulations. In Table I we have collected solvation frequencies obtained from a variety of previously reported simulations as well as from new calculations on model diatomic solvents. Mainly monatomic solutes with zero charge have been employed in the simulations listed. An exception are the "CHjCl" simulations of Carter and Hynes,' which involved a diatomic solute whose atoms were either uncharged (NP) or which carried a *le charge pair (IP). The solvation frequency listed for these simulations is that of the reaction potential difference between the twosolute sites. In thestockmayer solvents recently studied by Neria and Nitzangand Perera and Berkowitz,Io what was simulated was the nonequilibrium response (Si-&)) to a charge jump (0 le) rather than Cv(t) as in the other simulations. On the basis of results by these authors9and other^,^^^ we assume that the initial frequency of this response is identical to w, of Cv(t) of the initial, i.e. uncharged, solute. Finally, the entries labelled 'DD" represent new simulations of a series of 'dipolar dumbbell" solvents. These simulations involve an uncharged and stationary solute that interacts with 252 solvent molecules via site-site potentials of Lennard-Jones plus Coulomb type. The algorithm and boundary conditions used in these
-
The Journal of Physical Chemistry, Vol. 97, No. 1, 1993 15
Letters
v1
3
w
..VI,, 1 ,
0. 0.
8.
16.
24.
a
.o'
'
'
'
'
'
'
'
'
'
'
w
n -4 W
32.
Figure 1. Ratio of the solvation frequency usobserved in simulation to the free rotor frequency w1 plotted as a function of as.In all cases a6 was calculated from eq 13 neglecting the ( 1 - 1/60) term. Data are numbered according to Table I.
.o
.5
1.o
tu70 simulations are similar to those reported in ref 5. Solute LJ parameters are t/kg = 200 K, u = 3.5 A. The diatomic solvent consists of two sites separated by 2 A and having the same LJ parameters as the solute and 6+/& charges of magnitudes listed in Table I. The density and temperature employed were respectively 0.0134 A-3 and 298 K. In Figure 1 the values of ( W , / W ~ )obtained ~ from all of these simulations are plotted versus ascalculated on the basis of eq 13. Although the dielectric constants of most of the solvents have not been determined, with the exception of the O.le DD case, we expect co > 10 and have therefore neglected the (1 - 1/to) term in evaluating as.As can be seen from this figure and the tabulated data, agreement between the calculated and observed solvation frequencies is good. With only three exceptions, the predicted solvation frequencies usare within *lo% of the observed values. Two of the cases where poorer agreement is found involve small charged solutes in CH3CN. Further analysis indicates that eq 10 is still a good approximation in these cases but that eq 13 is less accurate due to the fact that such solutes induce significant changes in the solvent structure in their immediate neighborhood. For example, the solvent density is about 20% higher in the first solvation shells of these solutes than in the bulk. This higher dipole density means that the continuum estimate of (6P),is, slightly low and thus so too is the predicted us.However, even in these examples the predicted frequencies differ by no more than 20% from the observed values. The other instance where poor agreement is found is for the least polar DD solvent (O.le). Here the lack of agreement mainly reflects the fact that it is no longer appropriate to ignore the (1 - l/to) term in eq 13. If the value of €0 = 2.8 obtained from an MSA calculation12for a dipolar hard-sphere solvent with the same dipole density is assumed, one obtains a value of us= 2.46 ps-l, in excellent agreement with the simulated value (2.4 ps-I). On the basis of the above comparisons we therefore conclude that eq 13 yields quantitatively accurate predictions of the solvation frequency, at least in cases where the solvent structure is not greatly perturbed by the solute. In such instances the fast Gaussian part of the solvation dynamics can therefore be understood in terms of this very simple model connecting the free rotations of single solvent molecules and the solvation coordinate 6V. The Complete Response Thus far our discussion has concerned only that part of the solvation response due to free inertial dynamics of solvent molecules. However, evidence from previous simulations of a Brownian dipole lattice solvents suggest that the connection embodied in eq 13 is more general. The solvent in these latter simulations consisted of an array of point dipoles fixed to the sites
Figure 2. Results of simulations of uncharged solutes in Brownian dipole lattice solvents (see ref 8). In panel (a) the dashed curves are Cl(t) and the solid curves Cv(t). The arrows indicate increasing solvent polarity. In panel (b) the solid curves are the simulated Cv(t) functions and the dashed curves are thosecomputed from C l ( t )using eqs 15 and 10. Values ofcr,usedwere 1.30,1.58,2.26,3.65,5.30(seeref8).Curvesfordifferent polarities (increasing from top to bottom) have been vertically displaced for clarity.
of a simple cubic lattice. The equations governing their motion were solved in the overdamped or purely diffusive limit, where inertial motions are absent. In these simulations we observed that the 1/ e time for relaxation of the solvation tcf ("T,") in these systems could be directly related to the same fluctuation ratio appearing in eq 10. That is, under a wide range of conditions, we found that
where TO is the decay time of the single-particle reorientational tcf ( C l ( t ) ,eq 3) in the absence of dipole-dipole interactions? Since the same factor asapparently translates between singleparticle rotations and the solvation response in both the inertial and the diffusive regimes, it is reasonable to suspect that there might be an simple link between Cl(t)and Cv(t)approximately valid for all times. With no a priori justification we suggest the relation
The utility of this relationship can be readily judged on the basis of simulation data. Relevant comparisons are provided in Figures 2-4. In Figure 2a we show tcfs observed in a series of Brownian dipole lattice solvents. The dashed curves in panel (a) are the single-particle tcfs CI(t) and the solid curves are the corresponding solvation tcfs Cv(t). With increasing solvent polarity the rotational motions become more and more sluggish whereas solvation becomes faster and faster. As shown in Figure 2b, the relationship between these contrary trends is captured by eq 15 to within the uncertainties in these data. The agreement shown here is quite remarkable and goes far beyond the relationship between T~ and T Oreported previously. (We note that in making these comparisons we have used eq IO rather than the more approximate eq 13 since exact values of (6v2),, are available for these lattices. Equation 13 predicts a,values that are uniformly higher by IO%.) Figure 3 presents a test of eq I5 using previous simulation results on two realistic solvent models, acetonitrileSand water.6
-
16
Letters
The Journal of Physical Chemistry, Vol. 97,No. 1, 1993
CH, CN n
4 W
>
u
n -4 W 3
u
.o
.1
.2
.3
.4
.5
Time ( p s ) Figure 3. Tests of eq 15 using simulations of uncharged solutes in acetonitrile CH3CN5and water.6 In each panel the solid curve is Cv(r), the long dashed curve is C l ( t ) , and the short dashed curve is Cv(t) calculated according to eqs 13 and 15. Values of asused are those listed in Table I.
In both of these cases the free inertial component of the solvation dynamics (solid curves) dominates the response. On the time scale over which Cv(t) decays by -90% there is only a very modest decay of Cl(t) (long dashed curves), implying average reorientation angles of
