3 February 1995
CHEMICAL PHYSICS LETTERS ELSEVIER
Chemical PhysicsLetters 233 (1995) 138-144
Computer simulation of solvation dynamics in several model solvents Yi Lin, Charles D. Jonah Chemistry Division, Argonne National Laboratory, Argonne, IL 60439, USA
Received 28 July 1994; in final form 14 November 1994
The effect on solvation dynamics and energetics of the length and charge distribution of a model linear solvent has been determined. Solvent clustering around the ion is independent of the solvent molecule length if the solvent dipole terminus is attracted to the charged solute. In contrast, if they repel, the structure around the solute ion depends on the length of the solvent molecules and on the solute charge. The calculated solvation dynamics has two components and depends on solvent structure, The fraction of the fast relaxation component decreases with the length of the solvent molecules.
While the role of solvation in chemistry has long been recognized, only recently has it been possible to probe the time scale of the solvent relaxation process using direct measurements [ 1-11 ]. In such studies, one considers the change in the solvent polarization that arises after the sudden alteration of the charge distribution in a solute molecule. The discovery of nonexponential solvation dynamics, the role of solvent molecule structure, and solvation time scales that differ from the predictions of continuum theory have prompted new and more sophisticated theoretical approaches [ 12-16]. Both our theoretical and experimental studies have emphasized the importance of the molecularity of the solvent (non-continuum nature) and the structure of the solvent molecule [ 10]. The current studies were undertaken to extend our understanding of the experimental measurements of anion solvation dynamics in alcohol solutions. The details of these experiments have been reported elsewhere [ 10]. In brief, it was found that there were large 0009-2614/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSD10009-2614(94)01421-3
shifts in the optical absorption spectrum of the anion as a function of time for primary alcohols, and the magnitude of the shift was independent of the length of the alcohols. However for secondary alcohols, the spectral shift was much smaller. The final relaxed spectrum of the anion in primary alcohols is blue-shifted relative to the absorption of the anion in secondary alcohols. These results are the converse of the predictions of a continuum theory of solvation, where the solvation energy would depend on the dipole density. The dipole densities are similar for alcohols that have the same number of carbon atoms but are strongly dependent on the carbon atoms in the alcohol. We suggested that the observed spectral differences of the benzophenone anion in linear and branched alcohol solvents of the same carbon number arise from the difference in the packing of the OH bonds around the anion. Using Monte Carlo techniques, we showed that this suggestion was indeed reasonable. However, Monte Carlo calculations do not probe the kinetics of solvation - instead they only model the equilibrium solvation configurations [ 10]. We have extended the
Y. Lin, C D . J o n a h / C h e m i c a l Physics Letters 233 (1995) 138-144
study to calculate the solvation times of the system using a molecular dynamics model. Closely related to the present work is that of Maroncelli and co-workers, who simulated the solvation dynamics of ions in acetonitrile [ 17] and in a solvent that was modeled as a Brownian lattice [ 18]. These simulations provide a detailed physical picture of the molecular structure and dynamics for these solvent systems. Karim et al. studied the response of the solvent to an instantaneous change of the solute dipole [ 19]. All of these studies display one common characteristic - the relaxation processes display nonexponential or multi-exponential behavior. The goal of the present calculations was to understand how molecular lengths and charges would affect solvation. Because many different parameters were to be surveyed, it was necessary to use a simple linear model, which may be most closely described as a simplified nitrile-like molecule. Using a classical rigid body molecular dynamics simulation technique developed by Fincham [20,21], we computed the radial distribution function for several different size solvent molecules and the solvation time-response function,
(V( t) ) - (V(o~) ) c ( t ) = ( v ( 0 ) ) - (v(o~)) '
smallest system studied, the size was approximately 35 /~. The simulation employed cubic periodic boundary condition with the minimum image convention. Initially the solvent molecules were placed on a facecentered cubic lattice with a lattice parameter of 6.7/~. The systems were considered in both a constant energy and constant temperature ensemble. Details of the simulation strategies are given in Ref. . The classical equation of motion was solved using the leapfrog scheme developed by Fincham for linear rigid molecules [ 21 ]. Time steps of 2 fs were used in the integration. With this scheme, the fluctuations of the total energy and rotational temperature were less than 0.2% in most cases. For the initial equilibration simulation, occasional rescaling of the kinetic and rotational energy was required. All the simulations reported here were performed at temperatures of approximately 300 K. The solvent molecules $2, $3 and $4 were modeled by a linear rigid array of n atoms as shown in Fig. la. $3 might be described as a model of acetonitrile. The dipole was modeled as partial charges, q + and q - on adjacent atoms so that the solvent molecule was electrically neutral. The '$2' solvent molecule consists of two atoms, both charged, while '$3' and '$4' have the two charged atoms and one or two neutral atoms, respectively. The center-to-center distance between the atoms was 1.2 A. The solutes were treated as immobile, monatomic ions of varying charge and were always located in the center of the simulation box. The shortrange part of the intermolecular interaction was modeled by a Lennard-Jones potential between the atoms. The long-range Coulomb interactions between the dipole groups of partial charges of the solvent were accounted for using a spherical approximation to the Ewald summation as described by Adams and Dubey [ 22]. The interaction parameters for both solvent and solutes are given in Table 1 and the complete interaction energy U is o
where V(t) is the electrical potential at the solute charge site (reaction potential) at time t. We found that both the equilibrium solvation structure and the relaxation dynamics vary with the size and shape of the solvent molecules involved. The fast initial relaxation, which occurred computationally on the time scale of 0.1--0.2 ps, was smaller in simulations that used longer solvent molecules than in simulations where shorter solvent molecules were used.
2. Simulation method Molecular dynamics simulations were performed for the system of 248 solvent molecules and a single charged solute molecule in the center of the simulation box. The size of the simulation cell and the interaction parameters were selected to be consistent with our previous studies. Thus, meaningful comparison to the results of our previous Monte Carlo simulations could be made. The dimension of the simulation box depends on the length of the model solvent molecule. For the
E ,a E Aj
(2) where i and j are molecular indices and a and/3 are atomic indices. The symbol Ai (Aj) represents the set of all atoms of molecule i. Only pairwise solvent-solvent and solvent-solute interactions were considered.
Y. Lin, C.D. Jonah / Chemical Physics Letters 233 (1995) 138-144
4= 6= D i s t a n c e (A)
D i s t a n c e (A)
..... + 1
Fig. 1. The radial distribution function g(r) as a function of distance r (in/~). r is measured from the solute ion to the center of the solvent dipole. $2, $3, and $4, denote the solvent molecules of length 2, 3, and 4 atoms respectively. (a) For a solute anion, (b) for a solute cation, (c) for a fixed solvent molecule $3, with different charges on the solute ion.
Table 1 Characteristics of the different solute and solvent properties simulated. In these simulations tr (.~) and E / k are the size parameter and potential parameter for the Lennard-Jones interaction Solute properties
E~ k ( K )
No. of atoms
Q0 QI Q2 Q3
2.1 2.1 2.1 2.1
0 1.0 2.0 - 1
$2 $3 $4
3.5 3.5 3.5
190 190 190
2 3 4
Y. Lin, C.D. Jonah / Chemical Physics Letters 233 (1995) 138-144
3. Results and discussion The solvent response to a sudden change of the charge of a solute molecule can be characterized both by how fast the solvent molecules rearrange and by the equilibrated configurations of the solvent around the solute. Both the dynamics and equilibrated configurations can be calculated using molecular dynamics; however, the analysis of the calculations will be different depending on whether equilibrium configurations or solvent dynamics are to be determined. To determine the time dependence one must observe the dynamics resulting from charge creation with a given initial spatial distribution of solvent molecules. This must be done for many different initial distributions of solvent molecules and the dynamics from the different initial configurations must be averaged. Only a few initial solvent configurations need to be studied to determine equilibrium properties; however the system must equilibrate before one can begin averaging equilibrium properties over the resulting configurations.
3.1. Equilibrium configurations Both the orientation and the distance of the model solvent dipole from the solute ion have been probed. We have also carried out the calculations for both a solute anion and cation. Fig. 1 displays the radial distribution function (g (r)) for the three different solvent molecules: (a) around a solute negative ion; (b) a solute positive ion and; (c) for $3 with four different charges on the solute. The function g(r) is the number of solvent dipoles in a unit volume at a distance r from the solute ion relative to the number that would be in the same volume if the solution were homogeneous (for a homogeneous solution g(r) would be 1). The distances were measured from the center of the ion to the midpoint of the dipole. For convenience, the terminus of the solvent molecule was always positive while the solute species could be either positive, negative or zero. Fig. la shows that all three model solvent molecules have very similar radial distribution functions for a solute anion. Both the distance and height of the first peak of the radial distribution function are similar and suggest that the solvation structure near the solute ion is very similar. This prediction is consistent with the experimental data for the final spectrum of the benzophenone anion in alcohols, where the shifts of
the spectra are independent of the length of the solvent alcohol molecules [ 10]. For the positive ions (Fig. lb), there is a considerable change in the magnitude of the first peak of the radial distribution function even though there are only slight differences in the position of the peak. The smaller value of the peak of the radial distribution function shows that there are fewer dipoles near the solute anion. These differences can be ascribed to the steric hindrance between the large solvent molecules. The solvent molecules do not point directly towards the solute ion (see below) and will interfere with one another. Thus there will be fewer solvent molecules near the solute molecule. Fig. lc displays g(r) as a function of the solute ion charge for the '$3' solvent. The results show that a solute anion will have more $3 molecules in the first solvation shell (defined as the distance from the solute ion to the first minimum, near 4.4,, in the radial distribution function) than will a solute cation. Close examination of Fig. Ic suggests that the maximum of the radial distribution function is closer for a solute cation than for a solute anion. Presumably this occurs because of the repulsion between the positive end of the solvent molecule and the solute positive charge means that the solvent is approximately perpendicular to the line between the solute and the center of the dipole and the distance is measured to the center of the dipole and not the closest point of the molecule. Fig. 2 shows the average orientation of the solvent molecules in the first solvation shell, (cos( /x.r) ), where /~ is the unit vector of the dipole in the solvent molecule (pointing from the negative charge to the positive charge), and r is the unit vector from the solute to the center of the dipole of the solvent molecule. This average would be - 1 if all the nearest molecules pointed directly toward the solute ion, 0 if the molecules are perpendicular to the line from the solute ion and 1 if the positive end of the molecule points directly away from the solute ion. The Coulomb attraction causes the positive end of the solvent molecule to point towards the solute anion. The alignment around a solute cation is weaker and in the reverse direction. This result reflects the balance of energetic constrains- the advantage of close approach between the negative end of the solvent dipole (which is not at the end of the molecule) and the solute cation, the repulsion between the solute cation and the positive end of the solvent dipole and
Y. Lin, C.D. Jonah/Chemical Physics Letters 233 (1995) 138-144
-0.5 -20 -
Fig. 2. The average orientation of the solvent molecules in the first solvation shell ( c o s ( g - r ) ) , is plotted as a function of solute ion charge Q0 for several model solvents. (,It) represents the solvents with two atoms ($2), (IS]) represents $3, and (©) represents $4. The lines are drawn only to guide the eyes.
the short-range repulsion between the neutral atom (s) of the solvent molecule and the solute molecule. In the solvation experiments, one measures the shifts in the absorption and emission spectra of the charged solute molecule and not the alignment of the solvent molecules nor the energy of solvation. The shift in the electronic spectrum of the solute ion arises from the effect of the field created by the solvent dipoles on it. We can estimate this shift by calculating the average reaction potential (V), the electrical potential at the solute site. Fig. 3 displays the reaction potential as a function of solute charge for $2, $3, and $4. These results indicate that the potential is independent of the length of the solvent molecule for anions. This prediction is consistent with our experimental results [ 10], which show that the spectrum of the equilibrated solvated benzophenone anion is virtually identical for primary alcohols from n-propanol to n-decanol. From these results, we can infer that the local structure, rather than the long range interaction, dominate the solvation energetics for solute anions. The density of the solvents $2, $3 and $4 were kept the same in the calculation of Fig. 3. As a consequence, the dipole density of $4 will be less than the dipole density of $3 or $2. The potential at the solute anion is the same for all three solvents, showing that the potential does not depend on dipole density. However, as seen in Figs. 1
Fig. 3. The electric potential (V) at the solute ion site is plotted as a function of solute charge Q0 for solvent molecules of length 2, 3, and 4 atoms. ( O ) denotes $2, (IS]) denotes $3 and ( A ) denotes $4.
and 2, the local structure is the same for the different solvent lengths. Because the energetics remain the same for $2, $3 and $4 as do the local structure while the density of dipoles decreases, we believe that the effect of the local structure is dominant on the energetics at the solute molecule. We examined the reaction potential as a function of the solute charge to test the validity of the linear response assumption for this system. If the linear response assumption is valid, the potential energy at the solute would be proportional to the solute charge. Because of the linearity of (V) versus charge for a solute anion, the results in Fig. 3 show that the linear response theory may be adequate to predict the energy for a solute anion. However the major departure from linearity for central positive ions shows that this concept is not appropriate for a highly charged solute cations. Interestingly enough, for highly charged solute cations, the energies may be more consistent with a continuum model where the amount of solvation energy is dependent on the dipole density. The ratio of (V) values for $2, $3, and $4 for a central ion charge of + 2 is almost exactly 4 : 3 : 2; the ratio of the dipole densities. This may result from the lack of radial order around the solute cation.
3.2. Dynamics The response of the solvent to an introduction of the solute charge was computationally observed to study
Y. Lin, C.D. Jonah/Chemical Physics Letters 233 (1995) 138-144 1.0-
600 Time (fs)
Fig. 4. The solvation correlation function C(t) is plotted as a function of time for solute charges ( - 1, + 1 ) and solvent molecules of length $3, $4. The two-component relaxation dynamics are evident.
solvation dynamics. The value of C(t) (Eq. (1) ) was determined as a function of time, rather than using the linear response assumption and calculating C(t) from the equilibrium time-correlation function [ 17]. Before the calculation of the solvation dynamics, the solvent was equilibrated around a neutral solute molecule. Subsequently a charge was located on the solute molecule, and the kinetic energy, rotational energy, and (V(t)), the electrical potential due to the solvent molecules at the solute site, were calculated for approximately 2 ps. This process was repeated for 500 different initial configurations to study the solvation dynamics. The nonequilibrium solvation-time-response function C(t) was determined for the $3 and $4 molecule for a univalent cation and anion. The results are given in Fig. 4. Common to all the systems studied, the dynamics as measured by C(t) show two components, a fast response of 0.1-0.2 ps and a slower relaxation of 1 ps or greater. The fast component is nonlinear, with an initial slope of zero. As can be seen in Fig. 4, the solvation dynamics are slightly different for the relaxation around a cation and anion. The fraction of the fast relaxation is greater for the shorter solvent molecules than for the longer solvent molecules. The relaxation process around the solvent can be thought as taking place via three relaxation modes: inertial rotation, whereby the solvent molecule rotates towards the charged species without collision; rotational diffusion, where the solvent molecule rotates
towards the charged species but collides with the other solvent molecules; and translational diffusion, whereby the solvent molecule diffuses towards the solute molecule. If we assume that the dominant processes for solvation on this time scale are inertial rotation and diffusional rotation, the dependence of the duration and amount of the fast component on the size of the solvent molecule can be understood qualitatively. As the solvent molecule becomes longer, the probability of the molecule striking a neighboring solvent molecule during the rotation increases. This will change the motion from inertial rotational motion to a diffusive rotational motion. For a solute cation, the solvent molecules will not be as strongly aligned at equilibrium as they are for a solute anion. One might expect that less solvent motion would be necessary to solvate the cation solute and thus the inertial motion, which would remain about the same size, would account for a greater fraction of the relaxation. This is consistent with the results shown in Fig. 4. We have also examined the translational diffusion by calculating the mean-square displacement as a function of time for the $3 model. We found that diffusion was quite small for the model systems considered here. For this reason and the fact that the hypothesis of inertial rotation and rotational diffusion provided an explanation for our results as a function of length, we consider this possibility a less likely explanation.
4. Summary We have investigated the application of molecular dynamics simulation to study ion solvation in dipolar media. Our simple model is capable of qualitatively predicting the different spectral shifts that occur in the solvation of benzophenone anions in alcohols and in acetonitrile. For the solvation of an anion in an alcohollike solvent, the orientation, closest approach, and number of solvent molecules in the first solvation shell are the same, independent of the length of the solvent molecule. This leads to the potential at the solute being independent of the length of the solvent molecule, which means that the solvation energetics will be dominated by the first solvation shell. These results are in agreement with the experimental results, which show that the spectral shift of the benzophenone anion is independent of the length of the solvent molecule. Ear-
Y. Lin, C.D. Jonah/Chemical Physics Letters 233 (1995) 138-144
lier, we had studied solvation in clusters using M o n t e Carlo techniques and the importance o f local structure found in those p r e v i o u s studies has been reproduced using m o l e c u l a r d y n a m i c s [ 10]. The calculated d y n a m i c s display two components, as o b s e r v e d p r e v i o u s l y by M a r o n c e l l i and c o - w o r k e r s in simulations o f solvation o f a simple anion in acetonitrile [ 17]. W e h a v e s h o w n that the a m o u n t o f the fast relaxation c o m p o n e n t with a zero initial slope depends on the length o f the solvent m o l e c u l e . This is consistent with the fast c o m p o n e n t arising f r o m primarily inertial rotation o f the solvent m o l e c u l e w h i l e the slower term m a y be due to rotational diffusion ( a n d possibly translational diffusion) o f the solvent molecules. Further studies are needed on this point. W e h a v e b e g u n calculations on m o r e c o m p l i c a t e d solvent m o l e c u l e shapes to confirm the results obtained p r e v i o u s l y using our M o n t e C a r l o simulations for prim a r y and secondary alcohol m o l e c u l e s [ 10]. W e plan to study both solvation d y n a m i c s and energetics for these systems.
References [ 1] Solvation, Faraday Discussions Chem. Soc. 85 (1988).  G. Vander Zwan and J.T. Hynes, J. Phys. Chem. 89 (1985) 4181.
 P.G. Wolynes, J. Chem. Phys. 86 (1987) 5133.  I. Rips and J.J. Jortner, J. Chem. Phys. 87 (1987) 2090. [ 5 ] E.W. Castner, M. Maroncelli and G.R. Fleming, J. Chem. Phys. 86 (1987) 1090.  S.G. Su and J.D. Simon, J. Phys. Chem. 91 (1987) 2693.  V. Nagarajan, A.M. Brearley, T.J. Kang and P.F. Barbara, J. Chem. Phys. 86 (1987) 3183.  G.A. Kenney-Wallace and C.D. Jonah, J. Phys. Chem. 86 (1982) 2572.  A. Migus, Y. Ganduel, J.L. Martin and A. Antonetti, Phys. Rev. Letters 58 (1987) 1559. [ 10] Y. Lin and C.D. Jonah, J. Phys. Chem. 96 (1992) 10119; 97 (1993) 295. [ 11 ] F.H. Long, H. Lu and K.B. Eisenthal, J. Chem. Phys. 91 (1989) 4193. [ 12] M.A. Kahlow, T.J. Kang and P.F. Barbara, J. Chem. Phys. 88 (1988) 2372. [ 13] C.F. Chapman, R.S. Fee and M. Maroncelli, J. Phys. Chem. 94 (1990) 4929. [ 14] M. Maroncelli and G.R. Fleming, J. Chem. Phys. 89 (1988) 5044. [ 151 M. Maroncelli, V.P. Kumar and A. Papazyan, J. Phys. Chem. 97 (1993) 13. [ 16] B. Bagchi, Ann. Rev. Phys. Chem. 40 (1989) 115. [ 17] M. Maroncelli, J. Chem. Phys. 94 ( 1991 ) 2048.  A. Papazyan and M. Maroneelli, J. Chem. Phys. 95 (1991) 9219. [ 19] O.A. Karim, A.D.J. Hayment, M.J. Banet and J.D. Simon, J. Phys. Chem. 92 (1985) 3391. [ 20] M.P. Allen and D.J. Tildesley, Computer simulation of liquids (Oxford Univ. Press, Oxford, 1987).  D. Fincham, More on rotational motion of linear molecules, CCP5 Quarterly 12 (1984) 47.  D.J. Adams and G.S. Dubey, J. Comput. Phys. 72 (1987) 156.