An Investigation of Water Dynamics in Binary Mixtures of Water and

J. Phys. Chem. B 2008, 112, 283-293

283

An Investigation of Water Dynamics in Binary Mixtures of Water and Dimethyl Sulfoxide† Michael R. Harpham, Nancy E. Levinger,* and Branka M. Ladanyi* Department of Chemistry, Colorado State UniVersity, Fort Collins, Colorado 80523-1872 ReceiVed: June 26, 2007; In Final Form: October 3, 2007

The motion of water molecules in mixtures of water and d6-dimethyl sulfoxide (DMSO) has been explored through molecular dynamics (MD) simulations using the SPC/E water model (J. Chem. Phys. 1987, 91, 6269) and the P2 DMSO model (J. Chem. Phys. 1993, 98, 8160). We evaluate the self-intermediate scattering functions, FS(Q,t), which are related by a Fourier transform to the incoherent structure factors, S(Q,ω), measured in quasielastic neutron scattering (QNS) experiments. We compare our results to recent QNS experiments on these mixtures reported by Bordallo et al. (J. Chem. Phys. 2004, 121, 12457). In addition to comparing the MD data to the experimental signals, which correspond to a convolution of S(Q,ω) with a resolution function, we examine the rotational and translational components of FS(Q,t) and investigate to what extent simulation results for the single-molecule dynamics follow the dynamical models that are used in the analysis of the experimental data. We find that the agreement between the experimental signal and the MD data is quite good and that the portion of FS(Q,t) due to translational dynamics is well represented by the jump-diffusion model. The model parameters and their composition dependence are in reasonable agreement with experiment, exhibiting similar trends in water mobility with composition. Specifically, we find that water motion is less hindered in water-rich and water-poor mixtures than it is near equimolar composition. We find that the extent of coupling between rotational and translational motion contributing to FS(Q,t) increases as the equimolar composition of the mixture is approached. Thus, the decoupling approximation, which is used to extract information on rotational relaxation from QNS spectra at higher momentum transfer (Q) values, becomes less accurate than that in water-rich or DMSO-rich mixtures. We also find that rotational relaxation deviates quite strongly from the isotropic rotational diffusion model. We explore this issue further by investigating the behavior of orientational time correlations for different unit vectors and corresponding to Legendre polynomials of orders 1-4. We find that the rotational time correlations of water molecules behave in a way that is more consistent with the extended jump rotation model recently proposed by Laage and Hynes (Science 2006, 311, 832).

I. Introduction Mixtures of dimethyl sulfoxide (DMSO) and water are of considerable importance as chemical reaction media and in biomedical applications,1-4 which exploit the cryoprotection properties of these solutions.3,4 Given that interspecies hydrogen (H) bonds are formed between water hydrogen sites and the oxygen sites of DMSO, these mixtures are strongly associating and exhibit nonideal behavior, especially in their transport properties.3-12 Maximum deviations from ideality are reported to occur between 0.3 and 0.4 mole fraction DMSO (XDMSO). For example, the mixture has a higher viscosity than the pure components, with the peak occurring around XDMSO ) 0.33.5,6 These strong deviations from ideality have been hypothesized to arise from the formation of 2:1 water-to-DMSO complexes.8,11,13 A neutron diffraction study of water-DMSO mixtures14,15 found that water-water H bonding in mixtures was reduced relative to bulk water and that a large proportion of the water hydrogen atoms within the mixture were associated with the lone pairs on the DMSO oxygen. These data also suggest the presence of a high concentration of 1:2 DMSOwater complexes at the composition XDMSO ) 0.33.15 The †

Part of the “James T. (Casey) Hynes Festschrift”. * Corresponding author. E-mail: [email protected]; [email protected].

neutron diffraction data show that there exists local tetrahedral order of water molecules. A number of molecular dynamics (MD) simulations aimed at elucidating several aspects of the structure and dynamics of water-DMSO mixtures have been performed.16-29 Simulations performed by Luzar and Chandler17 show that the water-DMSO hydrogen (H) bonds live longer than water-water H-bonds, and that DMSO is usually a double H-bond acceptor that forms 2:1 aggregates of DMSO and water with almost tetrahedral symmetry. The MD simulation study by Borin and Skaf20 explored the local structures and H-bond distributions over the full composition range of water-DMSO mixtures. On the basis of their analysis of site-site pair distribution functions, they find that the water-DMSO hydrogen bonding is stronger than water-water hydrogen bonding. Furthermore, by adding DMSO to water, a hydrogen-bond-accepting water molecule is replaced by a DMSO molecule, which preserves the local tetrahedral order. They found that the 1:2 complexes of DMSO and water predominate over 2:1 complexes for water-rich mixtures (XDMSO < 0.5), and the opposite is true for water-poor mixtures (XDMSO > 0.5). They explored the self-diffusion coefficients and singleparticle rotational relaxation times of water and DMSO as functions of composition, and found that the water dynamics are slowest at the XDMSO ) 0.5 composition. Single-molecule dynamics of water and DMSO molecules in the mixtures have been studied by several spectroscopic

10.1021/jp074985j CCC: $40.75 © 2008 American Chemical Society Published on Web 12/05/2007

284 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Harpham et al.

Figure 1. SPC/E water (left) and P2 DMSO (right) models with structural parameters superimposed. Oxygen, carbon, hydrogen, and methyl groups are given in red, yellow, light blue, and dark blue, respectively.

TABLE 1: LJ and Coulomb Potential Parameters molecule

site

R/kB (K)

σR (Å)

qR (e)

DMSO

oxygen sulfur methyl oxygen hydrogen

35.988 119.96 147.93 78.177

2.8 3.4 3.8 3.166

-0.459 0.139 0.160 -0.8476 0.4238

water

TABLE 2: Composition and Molar Volumes of Water-DMSO Mixtures XDMSO

Nwater

NDMSO

molar volume (10-5 m3/mol)

0.124 0.25 0.33 0.50 0.66 0.75

438 375 335 250 165 125

62 125 165 250 335 375

2.52 3.12 3.50 4.37 5.37 5.81

methods, including Raman9 and NMR10,11,30 spectroscopy. A method that provides information about the translational and rotational mobility of molecules in the liquid phase is quasielastic neutron scattering (QNS). QNS experiments measure the incoherent dynamic structure factor, S(Q,ω), as a function of momentum transfer Q and energy transfer ω. Since the 1H isotope of hydrogen has a much larger incoherent cross section than D,31 selective deuteration can be used to obtain information on the dynamics on each of the mixture components. QNS has recently been applied to water-DMSO mixtures, first as a function of temperature at the eutectic composition13,32 (2:1 water-DMSO mole ratio) and then as a function of composition at room temperature.33 Specifically, in the latter study, Bordallo et al.33 reported results from QNS experiments on water/DMSO mixtures over a range of XDMSO from 0.1 to 0.8. In both studies, information about the rotational and translational dynamics was extracted from experimental data using the jump-translational diffusion and isotropic rotational diffusion models, which are well-established in their application to the analysis of QNS data from liquid samples.34,35 In the case of d6-DMSO/H2O mixtures, Bordallo et al.’s experiments indicated that the translational diffusion constant of water reaches a minimum of 0.3 × 10-5 cm2/s at XDMSO ) 0.35. This is about one order of magnitude smaller than the diffusion coefficient of bulk water.34 The rotational diffusion coefficient did not exhibit as much variation within this range, which seems to conflict with the MD simulation data20 as well as the NMR data.10 Estimates of the rotational correlation times of water obtained by NMR10 show a similar trend to the QNS data, but a larger variation with XDMSO. That is, both techniques indicate that XDMSO ) 0.33 shows the slowest rotational relaxation; however, the rotational relaxation time for the O-H vector reorientation varies over a larger range in the NMR experiments than it does in the QNS experiments. The rotational relaxation time of water O-H in the XDMSO ) 0.33 mixture,

Figure 2. Water total ISFs FS(Q,t) at Q ) 1.0 Å-1 (top) and Q ) 1.4 Å-1 (bottom), for water-DMSO mixtures of compositions XDMSO ) 0.124 (squares), 0.25 (filled circles), 0.33 (filled triangles), 0.5 (diamonds), 0.66 (open triangles), and 0.75 (open circles).

when measured by NMR,10 is approximately 4.2 times the value in neat liquid water, compared to 1.8 times when measured by QNS.33 MD simulation of the rotational relaxation contributions to S(Q,ω) can help elucidate the reasons for this discrepancy. In the work reported here, we perform a series of MD simulations on water-DMSO mixtures to gain information about structural and dynamical properties at a variety of compositions, with the aim of gaining a better understanding of the QNS experiments. We focus on the systems in which DMSO is deuterated and investigate the single-molecule dynamics of water. The primary quantity that we evaluate is the self-intermediate scattering function, which is related to the incoherent structure factor measured in QNS experiments. Specifically, S(Q,ω) is the frequency Fourier transform of the self-intermediate scattering function, FS(Q,t): S(Q,ω) )

1 2π

∫

∞

-∞

FS(Q,t)eiωtdt

(1)

where FS(Q,t) )

1 N

N

∑ 〈exp{iQ‚[r (0) - r (t)]}〉 j

j

(2)

j)1

N is the number of water hydrogen atoms, and rj(t) is the position of jth hydrogen atom. MD simulation provides a means of obtaining FS(Q,t) from molecular trajectory data. By evaluating the contributions to these self-intermediate scattering functions from translational and rotational dynamics, we can

Water Dynamics in Water-DMSO Mixtures

-1 Figure 3. Water translational ISFs FCM (top) and S (Q,t) at Q ) 1.0 Å Q ) 1.4 Å-1 (bottom), for water-DMSO mixtures of compositions XDMSO ) 0.124 (squares), 0.25 (filled circles), 0.33 (filled triangles), 0.5 (diamonds), 0.66 (open triangles), and 0.75 (open circles).

determine how these dynamics contribute to QNS and assess the range of validity of the dynamical models used in the analysis of experimental QNS data. In section II of this work, we present the model we use in our MD simulations, specify the time-correlation functions evaluated, and introduce the models that have been used in the interpretation of QNS data from water-DMSO mixtures. Our results and discussion are presented in section III. The MD results from FS(Q,t) are discussed in terms of rotational and translational relaxations of water molecules within the mixture, and comparisons between our data and the experimental QNS data obtained by Bordallo et al. are presented.33 The paper is concluded in section IV. II. Molecular Dynamics Simulations II.A. Interaction Model. Several combinations of interaction models have been compared in simulation studies of water/ DMSO mixtures.17,24,25 Of the model combinations in which the CH3 groups in DMSO are represented as single interaction sites, these studies indicate that the best results for the properties of mixtures at ambient conditions are obtained for the SPC/E water model36 and the P2 DMSO model.17 Given that in the present study we focus on the H2O/d6-DMSO system, a fully atomistic representation of DMSO is not needed to calculate the QNS spectrum of the mixture, thus we have chosen the SPC/E and P2 models for this study. A fully atomistic model25 of DMSO would be needed to simulate the QNS spectra corresponding to D2O/DMSO mixtures.

J. Phys. Chem. B, Vol. 112, No. 2, 2008 285

Figure 4. Water rotational ISFs FRS (Q,t) at Q ) 1.0 Å-1 (top) and Q ) 1.4 Å-1 (bottom), for water-DMSO mixtures of compositions XDMSO ) 0.124 (squares), 0.25 (filled circles), 0.33 (filled triangles), 0.5 (diamonds), 0.66 (open triangles), and 0.75 (open circles).

We have considered water/d6-DMSO mixtures with XDMSO ) 0.124, 0.25, 0.33, 0.50, 0.666, and 0.75 at ambient temperatures. Intermolecular interactions in the SPC/E and P2 models17,36 are described in terms of a Lennard-Jones (LJ) plus Coulomb form, with the interactions between sites R and β on different molecules given by uRβ(r) ) 4(Rβ)1/2

[(

) (

σ R + σβ 2r

12

-

)]

σR + σβ 2r

6

+

qRqβ 4π0r

(3)

LJ potential parameters, R and σR, and partial charges qR are summarized in Table 1. Figure 1 displays the bond angles and bond lengths for the P2 DMSO (1a) and SPC/E water (1b) molecules. The simulations were performed on NVE ensemble for systems of 500 molecules in a cubic box with periodic boundaries at an average temperature of 298 K. Box lengths were chosen to match the experimental densities7 at room temperature and ambient pressure. Compositions and molar volumes are listed in Table 2. LJ interactions were cut off at half the box length, while Ewald sums with conducting boundary conditions37 were used for Coulomb forces. Integration of equations of motion was accomplished using the velocityleapfrog algorithm for translations38 and the quaternion algorithm for rotations.39 Production trajectories were 0.5-1.5 ns, based on water content and expected mobility, with longer trajectories for lower water content and lower mobility compositions. Before the production trajectories, mixtures were equili-

286 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Harpham et al.

Figure 5. Water product and total ISFs FP(Q,t) and FS(Q,t) at Q ) 1.0 Å-1 (left) and Q ) 1.4 Å-1 (right), for water-DMSO mixtures of compositions XDMSO ) 0.25 (circles) and 0.33 (triangles) in the top graphs, and 0.66 (triangles) and 0.75 (circles) in the bottom graphs. Filled symbols represent total ISFs, whereas open symbols represent the product approximation to the ISFs.

brated starting from a face-centered cubic (fcc) structure, with molecules of different species randomly placed at lattice sites, with random orientations, at ∼15% higher molar volume, decreasing in steps to the experimental molar volume5 over 100 ps, followed by ∼100 ps of further simulation. In some cases, such as XDMSO ) 0.50, additional heating and cooling cycles were necessary to successfully equilibrate the mixture. II.B. Evaluation of the QNS Intermediate Scattering Function. The intermediate scattering function (ISF) FS(Q,t) can be obtained from MD simulation through eq 2. We focus our calculations on the dynamics of water molecules and the Q-range of the experiment on H2O/d6-DMSO mixtures performed by Bordallo et al.33 The time scale of a QNS experiment is determined by the resolution of the instrument. Because experimentally accessible time scales are typically on the order of tens of picoseconds, we calculate time correlations over similar time intervals. The QNS experiments of Bordallo et al. were performed on the QuasiElastic Neutron Spectrometer (QENS) instrument at Argonne National Laboratory, which has an instrument resolution of 90 µeV, corresponding to a temporal window of ∼50 ps, and a Q-range of 0.38 Å-1 < Q < 2.63 Å-1. In our simulations, we calculate the Q values in the same range at increments that are determined by the length (b) of the cubic box used in the simulations; that is, in each Cartesian direction, the Q values are integer multiples of 2π/b. In our simulations, we use a rigid model for water. We can thus predict the contributions to QNS from rotational and translational, but not from vibrational, motions of water molecules. This is not a serious drawback, given that vibrations in water occur at high frequencies that exceed the limits of detection in QNS.35 As a result, they contribute to the water QNS spectrum S(Q,ω) only through a frequency-independent multiplicative Debye-Waller factor.35

The location of a hydrogen atom in a water molecule can be expressed as r ) rCM + d

(4)

where d is the distance of the hydrogen atom from the molecular center-of-mass, rCM. In the rigid SPC/E water model, the molecular geometry is fixed and therefore d changes only through molecular rotation. The translational component of the water motion is tracked through rCM. The total ISF for a hydrogen atom in a water molecule is then given by FS(Q,t) ) 〈exp{iQ‚[rCM(0) - rCM(t)]} exp{iQ‚[d(0) d(t)]}〉 (5) In real systems, molecular rotational and translational motions are coupled. However, QNS experiments are frequently analyzed by neglecting rotation-translation coupling, and eq 5 would be factored into separate averages over rotational and translational coordinates. To test this approximation for the water-DMSO mixtures and to make contact with the experiments of Bordallo et al.,33 we compare the MD results for FS(Q,t) in the presence and in the absence of coupling. Assuming independence of rotational and translational dynamics in eq 5 results in the product approximation R FS(Q,t) = FP(Q,t) ) FCM S (Q,t)FS (Q,t)

(6)

where the translational (CM) and rotational (R) ISFs are given by FCM S (Q,t) ) 〈exp{iQ‚[rCM(0) - rCM(t)]}〉

(7)

Water Dynamics in Water-DMSO Mixtures

J. Phys. Chem. B, Vol. 112, No. 2, 2008 287 this dynamical model in their data analysis. According to this model, FCM S (Q,t) and Cl(t) decay exponentially, specifically -Γ(Q)t FCM S (Q,t) = e

(11)

Cl(t) = exp[-l(l + 1)DRt] ≡ exp[-t/τl]

(12)

and

where DR is the rotational diffusion constant, and τl is the rotational relaxation time for a Legendre polynomial of order l. As can be seen from eq 12, τl ) [l(l + 1)DR]-1 within this model. The decay constant of the translational ISF in the jumpdiffusion model41 is given by Γ(Q) )

DTQ2 1 + DTQ2τ0

(13)

where DT is the translational diffusion constant, and τ0 is the time between jumps. These two quantities are related through the jump length L by DT )

Figure 6. Orientational TCFs C1(t) (top) and C1,µ(t) (bottom) for water in water-DMSO mixtures of compositions XDMSO ) 0.33 (filled triangles), 0.5 (diamonds), and 0.66 (open triangles).

S(Q,ω) ) [j0(Qd)]2

1

∑ l)1

While the rotational ISF, FRS (Q,t), can be evaluated from MD trajectory data directly using eq 8, one can also express it in terms of the length d of the vector d and its orientation, dˆ ) d/d, to compare with models for rotational relaxation. This is done by using the Rayleigh expansion of the exponential40 ∞

FRS (Q,t) ) [j0(Qd)]2 +

∑ (2l + 1)[j (Qd)] C (t) 2

l

l

(9)

l)1

where jl is the spherical Bessel function of order l and Cl(t) is the time correlation for the Legendre polynomial Pl of order l of the change in direction of unit vector dˆ : Cl(t) ) 〈Pl[dˆ (0)‚dˆ (t)]〉

(10)

To extract information on rotational and translational dynamics from QNS data, it is convenient to use dynamical models that relate the results to familiar parameters such as the rotational and translational diffusion coefficients. Given that molecular motion in liquids is diffusive on intermediate and long time scales, a dynamical model appropriate for diffusive relaxation is usually applied to liquid-state QNS data. For bulk liquid water, a combination of the jump diffusion model for translation and the isotropic diffusion model for rotation fits experimental data reasonably well.35 In the present case, Bordallo et al.33 used

Γ(Q)

+

π [Γ(Q)]2 + ω2

∞

(8)

(14)

By comparing the ISFs from MD to the functional forms corresponding to the translational jump diffusion-isotropic rotational diffusion models, we can determine to what extent water dynamics in its mixtures with DMSO follow the behavior predicted by these models. Within the models for translational and rotational relaxation discussed above, S(Q,ω) is a sum of Lorentzians:

and FRS (Q,t) ) 〈exp{iQ‚[d(0) - d(t)]}〉

〈L2〉 6τ0

(2l + 1)

[jl(Qd)]2

Γ(Q) + l(l + 1)DR

π

[Γ(Q) + l(l + 1)DR]2 + ω2

(15)

Model parameters are obtained by fitting the experimental results to this functional form. The QNS spectrum also depends on molecular vibrations, which give rise to a multiplicative DebyeWaller factor, exp(-〈u2〉Q2/3). This factor, which depends on the mean-squared vibrational displacements, 〈u2〉, modulates the intensity of S(Q,ω), but does not affect its frequency dependence. Equation 15 shows that the low-Q spectrum is dominated by the translational contribution, given that [jl(Qd)]2 ∝ (Qd)2l at small Qd. At higher Q, the DR-dependent terms become important and increase the width and intensity of the spectral wings. Because the rotations contribute less to signals at lower Q than at higher Q, the assumption that translation and rotation are decoupled should show less effect at low Q than at high Q. Although eq 15 includes an infinite sum of Lorentzians, in practice, no more than four terms in the sum over l are needed to obtain convergence in the accessible range of Q values. III. Results and Discussion III.A. Computer Simulation Results. We explore the dynamics of water in these mixtures primarily through evaluation of the self-intermediate scattering functions, FS(Q,t). Figure 2 displays the total ISF, FS(Q,t), for the water in the mixtures at an experimentally relevant range of momentum transfer, Q.33

288 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Harpham et al.

Figure 7. First 10.0 ps (left) and 0.5 ps (right) of orientational TCF Cl(t) for l ) 1-4 for water in water-DMSO mixtures of compositions XDMSO ) 0.25 (top), 0.5 (center), and 0.75 (bottom).

We observe the slowest water dynamics as the mole fraction of DMSO approaches XDMSO ) 0.5. The same trend is observed for the center-of-mass ISF, FCM S (Q,t) (Figure 3), and the rotational ISF, FRS (Q,t) (Figure 4). Closer inspection of Figure 3 shows that translational motion appears nonexponential at short times ( 0.124 that

Water Dynamics in Water-DMSO Mixtures

J. Phys. Chem. B, Vol. 112, No. 2, 2008 289 diffusion model by investigating the behavior of the rotational TCF Cl(t), given by eq 10, for different l values. The isotropic rotational diffusion model assumes that water is a spherical rotor, that is, that its rotational diffusion tensor reduces to a single scalar constant, DR, and that molecular reorientation takes place through random, small-angle rotations. We test both of these assumptions. Because water molecules are not spherical, we anticipate that rotational relaxation might not be isotropic. Anisotropy in rotational relaxation can be tested by comparing the relaxation properties of different unit vectors embedded in the molecule. In this case, we compare the relaxation rate of the vector d to the unit vector uˆ along the molecular dipole. Comparison of Cl(t) with Cl,µ (t) ) 〈Pl[uˆ (0)‚uˆ (t)]〉

Figure 8. Results from fits of the long-time portion of the orientational TCF Cl(t) for l ) 1-4 to a single-exponential decay. The water orientational relaxation rates, 1/τl, are plotted vs l for water-DMSO mixtures of compositions XDMSO ) 0.124 (squares), 0.25 (filled circles), 0.33 (filled triangles), 0.5 (diamonds), 0.66 (open triangles), and 0.75 (open circles) in the top panel. The bottom panel contains τ1/τl plotted against l for the water-DMSO mixtures. Also included are the predictions of the isotropic rotational diffusion model (×) and of the Ivanov model with θ0 ) 60° (f).

could not be resolved through the QNS experiments of Bordallo et al.,33 in which the instrument resolution was on the order of 50 ps. We can use our MD data to determine the extent of rotationtranslation coupling present in the water ISFs. Figure 5 compares the full ISF, eq 5, to the product approximation, FP(Q,t), given by eq 6, for several compositions and Q values. The results displayed in Figure 5 reveal that the decoupling approximation is not consistently accurate within the Q range relevant to the experiment, which is 0.38 Å-1 < Q < 2.63 Å-1. From eq 15, we expect that the effects of coupling between translational and rotational motion on FS(Q,t) will increase with Q, and this result is seen in the MD data. In addition, the coupling between translation and rotation increases as the XDMSO ) 0.5 mixture is approached from either water-poor or water-rich compositions. This is consistent with the expectation that the coupling between water translational and rotational motions should increase as the hydrogen-bonding network becomes stronger and rotational motion becomes more hindered. Thus our results indicate that the decoupling approximation typically used in the analysis of experimental QNS data is not accurate for water-DMSO mixtures at high-Q and is least valid for the XDMSO ) 0.5 mixture. We further explore the rotational dynamics of water in these mixtures and evaluate the effectiveness of the isotropic rotational

(17)

allows us to determine the extent of rotational anisotropy of the system. In Figure 6 we display a comparison of Cl(t) and Cl,µ(t) for l ) 1 at several values of XDMSO. The figure shows that the relaxation rates of the water bond and dipole vectors are similar in water-rich mixtures (XDMSO e 0.33), but that, at higher XDMSO, their relaxation rates exhibit larger differences, with Cl,µ(t) showing faster decay than Cl(t). This indicates that rotational relaxation of water in mixtures with DMSO is not isotropic. The fact that relaxation in mixtures with XDMSO g 0.5 can be understood in terms of different effects on the bond and dipole vectors of the presence of water-DMSO hydrogen bonds, which are stronger than water-water hydrogen bonds.14,17,19,20 For water molecules bound to DMSO, the bond vector should relax more slowly than the dipole, given that its reorientation necessarily involves breaking these stronger hydrogen bonds. Since the water dipole vector can rotate through a relatively large angle before the hydrogen bond needs to break, the dipole vector relaxation is faster than O-H bond (actually, center-ofmass-H) vector relaxation. Anisotropy in water molecule rotational relaxation in neat liquid water has been observed previously in NMR experiments46 and MD simulations.46-48 Ludwig et al.12 have also reported NMR measurements of anisotropic water reorientation in water/DMSO mixtures for XDMSO ) 0.33 and found that the ratio of τ2 values for O-H and out-of-plane vectors increased with decreasing temperature. Figure 6 also shows that, as XDMSO approaches equimolar concentration, water rotational dynamics become more hindered, suggesting that interactions between water and DMSO molecules have the largest effect on the H-bond network strength and connectivity in the XDMSO ) 0.5 mixture. This trend agrees with the MD simulation studies performed by Borin and Skaf19,20 in which the results for Cl,µ(t) (l ) 1,2) at several values of XDMSO were reported. To further explore the rotational relaxation contributing to QNS, we investigate the l-dependence of the Cl(t) TCFs. Figure 7 displays plots of ln Cl(t)/[l(l + 1)] versus t at three values of XDMSO allowing results for l ) 1-4 to be compared. According to the isotropic rotational diffusion model, Cl(t) is given by eq 12, so a plot of ln Cl(t)/[l(l + 1)] as a function of t should be linear and independent of l. However, we observe that rotational relaxation becomes slower as l increases. This behavior is indicative of a high-torque liquid, and is expected for water.49 The right panels of Figure 7 demonstrate the behavior of ln Cl(t)/[l(l + 1)] at short times, 0 e t e 0.5 ps, which shows that this portion of Cl(t) is dominated by inertial and librational dynamics. As XDMSO increases, Cl(t) exhibits a larger number

290 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Harpham et al.

Figure 9. Comparisons between QNS spectra obtained from our MD simulations (lines) and QNS experiments (points) performed by Bordallo et al.33 at Q ) 0.7 (left), 1.0 (middle), and 2.6 Å-1 (right) for mixtures XDMSO ) 0.33 (top), 0.5 (middle), and 0.75 (bottom).

of librational oscillations, indicating that the librational motion of water becomes less damped as the number of water-DMSO hydrogen bonds increases. The oscillation frequency increases slightly with XDMSO, another manifestation of the fact that water-DMSO hydrogen bonds are stronger than water-water hydrogen bonds. Given that the l-dependence of the decay rate of ln Cl(t)/[l(l + 1)] is quite strong, a different model than rotational diffusion seems appropriate. For t > 0.5 ps, the decay is approximately exponential, decreasing with increasing l. A simple model that leads to this type of behavior is the rotational jump-diffusion model proposed by Ivanov.50 In this model, the molecule rotates in jumps rather than continuously. The molecule has a particular orientation during a residence time τ0, and then rotates instantaneously through an angle θ0. The relaxation time describing the time exponential decay of Cl(t) within the Ivanov model is given by

(

τjump ) τ0 1 l

sin[(l + 1/2)θ0] (2l + 1)[sin(θ0/2)]

)

-1

(18)

When θ0 becomes small enough so that the small-angle expansion of sin x (sin x = x) is valid, the Ivanov model reduces to the rotational diffusion model. Given that “jumps” through angles that are consistent with the hydrogen bond network geometry are expected in the present case, θ0 = 60° would be a reasonable estimate of the θ0 value.45 We have fit our Cl(t) data to single-exponential decay (Cl(t) ∝ exp(-t/τl)) for t g 1 ps. The results of this fit are displayed in Table 3 and Figure 8. The composition dependence of the τl values reflects the effects of the mixture hydrogen bond strength and connectivity

on water rotational relaxation. As can be seen from Table 3 and Figure 8, the slowest relaxation corresponds to water in the equimolar mixture, as anticipated. The values of the ratios τ1/τl, also shown in Table 3 and Figure 8, provide a simple test of the rotational relaxation model. In the rotational diffusion model, eq 12, the value is l(l + 1)/2, whereas, in the Ivanov model, this ratio depends also on the jump angle θ0:

τ1/τl )

τjump 1

) jump

τl

( 21)θ ]

[

3 (2l + 1) sin(θ0/2) - sin l +

0

(2l + 1)[3 sin(θ0/2) - sin(3θ0/2)]

(19)

For example, choosing θ0 ) 60° gives τjump 1 τjump l

)3-

6 sin[(2l + 1)π/6] 2l + 1

(20)

which leads to values of 2.40, 3.43, and 3.67 for l ) 2, 3, and 4, respectively. The corresponding values of the τ1/τl ratio are 3, 6, and 10 for the rotational diffusion model. We see from the results shown in the last two columns of Table 3 and in the bottom panel of Figure 8 that τ1/τl ratios fall significantly below the rotational diffusion model predictions at all mixture compositions. One set of results, for XDMSO ) 0.66, agrees quite well with the Ivanov model for θ0 ) 60°, as can be seen from Figure 8. However, our attempts to fit our results to this model at other compositions were less successful. We expect that these difficulties were due in part to a deficiency of this model that was recently pointed out by Laage and Hynes.45 Specifically, the model assumption that the molecule

Water Dynamics in Water-DMSO Mixtures

J. Phys. Chem. B, Vol. 112, No. 2, 2008 291

Figure 11. Water translational diffusion constants predicted from our fits of simulation mean-squared displacements (squares) and ISFs (circles), and QNS experiments (triangles) performed by Bordallo et al.33 Lines are drawn to guide the eye.

Figure 10. Example fits of the exponential decay constants Γ(Q) of the MD translational ISFs to the jump-diffusion model. The points are the results of fits of the decays of FCM S (Q,t) at specific Q values, while the lines are fits to eq 13.

remains stationary between jumps is not reasonable for hydrogenbonded liquids, given the fact that the hydrogen bond network is dynamic. Laage and Hynes45 proposed an extended jump model for reorientation of water in which hydrogen bond orientational relaxation is taken into account. This is done by including the relaxation rate of the oxygen-oxygen vector that connects the molecule of interest to its hydrogen-bonded partner. Thus the relaxation time for dˆ is given by 1 1 1 ) + τl τjump τOO l

(21)

l

where τOO is the oxygen-oxygen vector orientational relaxl ation time. This model has been successfully applied to SPC/E water at room temperature,45 using as input τOO calculated from the l orientational TCFs of the O-O unit vector determined from MD simulation. Application to water-DMSO would require a modification of the model to take into account the fact that different O-O vector relaxation times and different τjump l values are expected for water-water and water-DMSO pairs,

given the larger steric bulk of DMSO molecules and different strengths of water-water and water-DMSO hydrogen bonds. A similar situation has been considered recently by Laage and Hynes in their study of water orientational relaxation in NaCl aqueous solutions.51 Application of an extended jump rotational diffusion model to aqueous liquid mixtures such as water-DMSO would be an interesting task for future work and we expect that the resulting model would provide an improved representation of the ldependence of Cl(t). It should be noted that the presence of more than one type of hydrogen bond is expected to lead to deviations from single-exponential decay of the resulting orientational TCFs51 and that differences between the bulk XDMSO and the local DMSO mole fraction27 in the first coordination shell of a water molecule would help determine the relative contributions of dynamics of the two types of hydrogen bonds to water orientational relaxation. III.B. Comparisons to Experiment. Direct comparison of our simulations to experimental results allows us to determine whether the interaction models for water and DMSO adequately represent the translational and rotational dynamics observed in QNS. Convoluting our MD simulation results for S(Q,ω) with the approximately Gaussian experimental QNS resolution function Res(Q,ω), obtained by Bordallo et al.,33 permits direct comparison of simulation with experiment. Representative plots of S(Q,ω) X Res(Q,ω) at different Q and XDMSO values are shown in Figure 9. The overall agreement between MD simulation and experiment is very good, although the level deteriorates somewhat as the concentration of DMSO increases. With increasing XDMSO, the spectral wings predicted by MD have somewhat lower intensity than that observed experimentally. This indicates that water mobility in these mixtures is underestimated by the potential models used in the simulation. Since the level of agreement decreases with increasing DMSO concentration, the likely reason is that the water-DMSO potential is more attractive than it should be. The Lorentzian widths of low-Q data obtained experimentally agree well with the jump-translational diffusion model.33 We apply this model to our MD results for FCM S (Q,t). Because the jump-translational diffusion model assumes that single time scale for the translational motion, and the MD simulation shows more than one time component, as can be seen in Figure 3, we needed

292 J. Phys. Chem. B, Vol. 112, No. 2, 2008

Harpham et al. Figure 12 shows a comparison between the experimental and MD water translational residence times as a function of XDMSO. These have very similar composition dependence, and the values are fairly close. Although the translational motion of water is effectively captured by the jump diffusion model, our results indicate that the isotropic rotation model is not applicable for water dynamics in water-DMSO mixtures, as discussed in section III.A. Given that the MD results for Cl(t) are inconsistent with the ldependence of τl predicted by eq 12, the rotational relaxation parameters obtained from QNS experiments33 cannot be compared to the l-dependent relaxation times obtained from MD. It might be interesting to reanalyze the experimental QNS data for water in water-DMSO mixtures using the extended rotational jump-diffusion model, which appears to be more consistent with the behavior of Cl(t) found in our MD simulations.

Figure 12. Water residence time as predicted by our simulations (circles) and QNS experiments (triangles) performed by Bordallo et al.33

to consider how to fit the simulated data. We have excluded the short-time data (t < 5 ps) from the fit because dynamics are ballistic at very short times and sensitive to local structure for a subsequent short time scale. FCM S (Q,t) data in the window approximately from 5 to 80 ps are fit to a single-exponential decay with a relaxation time 1/Γ, where Γ would be the Lorentzian line width in the energy domain of a QNS scan. Figure 10 shows that, within these assumptions, the jumptranslational diffusion model (eq 11) fits the data reasonably well. A Lorentzian line width of less than -2 × 10-2 ps-1 would be narrower than the instrument resolution function in the experiments of Bordallo et al.33 Translational mobility can also be quantified in terms of mean-squared displacements of the molecules in the mixture. At long times, mean-squared displacements can be related to the translational diffusion coefficient, Deff, by the Einstein relationship: 〈|rCM(t) - rCM(0)|2〉 Deff ) lim tf∞ 6t

(22)

From the slope obtained by a linear fit to the mean-squared displacements in a time window from 5 to 15 ps, we extract the translational diffusion coefficient Deff, so named to preserve the nomenclature used by Borin and Skaf20 in their MD simulations of the water-DMSO system. Deff can be compared directly to DT, obtained by the jump-diffusion model. The resulting translational diffusion constants predicted from the mean-squared displacements and jump-translation fit of the MD data, are plotted in Figure 11 along with the experimental results from Bordallo et al.33 The translational diffusion constant predicted by the jump-diffusion analysis of the MD data agrees with that produced by a linear fit of the water mean-squared displacements, lending further support to the jump-translation model. The QNS33 and MD data for DT show similar composition dependence, but the values of the experimental diffusion coefficients are somewhat higher than the DT values obtained from MD. The difference is largest near equimolar composition, suggesting that the strength of the water-DMSO hydrogen bonds might be overestimated in the models used in our simulation. A similar observation was made by Borin and Skaf20 who compared their MD data to diffusion coefficients measured by NMR.30

IV. Summary and Conclusions This paper presents the single-particle dynamics of water in water/DMSO mixtures calculated from a set of MD simulations in the composition range relevant for comparisons to QNS experiments performed by Bordallo et al.33 We used the SPC/E water and the P2 DMSO models to carry out the simulation. To make contact with QNS experiments, the primary focus of our study was the single-molecule ISF, FS(Q,t). This quantity was compared to experiment and analyzed in terms of separate contributions from water rotational and translational dynamics as well as the extent of rotation-translation coupling. We summarize here our main findings. We find that the translational and rotational mobilities of water molecules are slowest in the equimolar water/DMSO, which agrees with other simulations of this system19-21 using the same model potentials. The combination of SPC/E water and P2 DMSO models may overestimate the water-DMSO attractive interactions,20 resulting in a growing discrepancy between diffusivities predicted by MD and extracted from experimental data as the composition approaches the XDMSO ) 0.5 mixture. To compare our MD simulations with the QNS experiments of Bordallo et al.,33 we calculated S(Q,ω) by Fourier transforming FS(Q,t) and constructed the convolution S(Q,ω) X Res(Q,ω) by using a fit to the experimental resolution function, Res(Q,ω). This comparison shows overall very good agreement with QNS data, although the MD dynamic structure factors become somewhat narrower than the experimental ones in DMSO-rich mixtures, indicating slower than experimental relaxation. A similar level of agreement with experiment extends to the predicted translational dynamics within the jump-translational model. We compared FS(Q,t) to its approximation as a product of R translational and rotational ISFs FCM S (Q,t) and FS (Q,t) to determine the range of validity of the decoupling approximation for rotations and translations, which was used in the analysis of QNS data by Bordallo et al.33 We found that the extent of coupling between rotation and translation varies with XDMSO and Q. At low-Q, where rotational effects are negligible, the approximation is reasonable over the full composition range. However, at high-Q, the approximation breaks down. The extent of hydrogen bonding also affects the coupling between translation and rotation. Stronger hydrogen bonding as seen at XDMSO ) 0.33, 0.50, and 0.66 leads to stronger coupling between water rotational motion and translational motion. This means that the decoupling approximation is least applicable at XDMSO ) 0.5,

Water Dynamics in Water-DMSO Mixtures and most reasonable in low-Q and water-rich or water-poor samples. Experimental results33 for translational motion are based on the fitting of low-Q data, where the decoupling approximation has negligible effect because the rotational contribution is relatively small, and is thus quite reasonable. However, rotational motion is analyzed at high-Q, where the decoupling approximation has the greatest effect, and thus is less reliable. An additional complication in extracting information on rotational relaxation from QNS comes from the fact that FSR(Q,t) is a weighted sum of orientational time correlations, Cl(t). Thus a model for the dependence of the relaxation rate on the order of l of the Legendre polynomial for the reorientation of the center-of-mass-H unit vector dˆ is needed. Our simulations show that the l-dependence of Cl(t) at all mixture compositions is significantly weaker than predicted by the isotropic rotational diffusion model used in the analysis of experimental data, but usually stronger than that for the Ivanov jump-rotation model.50 This type of l-dependence is consistent with the predictions of the extended jump model, recently proposed for neat liquid water.45 By comparing the rotational relaxation rates of dˆ and the molecular dipole, we found that water reorientation is not isotropic and that the extent of anisotropy increases with XDMSO. It would be interesting to find out if this trend can be experimentally confirmed, for example by using NMR, which makes it possible to measure relaxation times of more than one unit vector embedded in the molecule.12,46 Acknowledgment. This article is based upon work supported in part by the Department of Energy Grant No. DE-FG0302ER15376 to N.E.L. and B.M.L. and by NSF grants CHE 0608640 (B.M.L.) and CHE 0415260 (N.E.L.). We acknowledge the assistance of Mr. Bharat Viswamani of the Department of Electrical and Computer Engineering at Colorado State University who performed preliminary simulation studies of waterDMSO mixtures. References and Notes (1) Biological Actions of Dimethyl Sulfoxide; Jacob, S. W., Herschler, R., Eds.; New York Academy of Science: New York, 1975; Vol. 243, p 508. (2) de la Torre, J. C. Ann. N. Y. Acad. Sci. 1983, 411, 293. (3) Yu, Z. W.; Quinn, P. J. Biosci. Rep. 1994, 14, 259. (4) Hernandez-Perni, G.; Leuenberger, H. Eur. J. Pharm. Biopharm. 2005, 61, 201. (5) Cowie, J. M.; Toporowski, P. M. Can. J. Chem. 1961, 39, 2240. (6) Schickman, S. A.; Amey, R. L. J. Phys. Chem. 1971, 75, 98. (7) Aminabhavi, T. M.; Gopalakrishna, B. J. Chem. Eng. Data 1995, 40, 856. (8) Kaatze, U.; Pottel, R.; Schafer, M. J. Phys. Chem. 1989, 93, 5623. (9) Baker, E. S.; Jonas, J. J. Phys. Chem. 1985, 89, 1730. (10) Gordalla, B. C.; Zeidler, M. D. Mol. Phys. 1986, 59, 817.

J. Phys. Chem. B, Vol. 112, No. 2, 2008 293 (11) Gordalla, B. C.; Zeidler, M. D. Mol. Phys. 1991, 74, 975. (12) Ludwig, R.; Farrar, T. C.; Zeidler, M. D. J. Phys. Chem. 1994, 98, 6684. (13) Cabral, J. T.; Luzar, A.; Teixeira, J.; Bellissent-Funel, M. C. J. Chem. Phys. 2000, 113, 8736. (14) Soper, A. K.; Luzar, A. J. Chem. Phys. 1992, 97, 1320. (15) Soper, A. K.; Luzar, A. J. Phys. Chem. 1996, 100, 1357. (16) Vaisman, I. I.; Berkowitz, M. L. J. Am. Chem. Soc. 1992, 114, 7889. (17) Luzar, A.; Chandler, D. J. Chem. Phys. 1993, 98, 8160. (18) Liu, H.; Mueller-Plathe, F.; van Gunsteren, W. F. J. Am. Chem. Soc. 1995, 117, 4363. (19) Borin, I. A.; Skaf, M. S. Chem. Phys. Lett. 1998, 296, 125. (20) Borin, I. A.; Skaf, M. S. J. Chem. Phys. 1999, 110, 6412. (21) Skaf, M. S. J. Phys. Chem. A 1999, 103, 10719. (22) Benjamin, I. J. Chem. Phys. 1999, 110, 8070. (23) Kirchner, B.; Searles, D. J.; Dyson, A. J.; Vogt, P. S.; Huber, H. J. Am. Chem. Soc. 2000, 122, 5379. (24) Vishnyakov, A.; Lyubartsev, A. P.; Laaksonen, A. J. Phys. Chem. A 2001, 105, 1702. (25) Strader, M. L.; Feller, S. E. J. Phys. Chem. A 2002, 106, 1074. (26) Chalaris, M.; Samios, J. J. Mol. Liq. 2002, 98-99, 399. (27) Mu¨ller, M. G.; Hardy, E. H.; Vogt, P. S.; Bratschi, C.; Kirchner, B.; Huber, H.; Searles, D. J. J. Am. Chem. Soc. 2004, 126, 4704. (28) Mancera, R. L.; Chalaris, M.; Samios, J. J. Mol. Liq. 2004, 110, 147. (29) Nieto-Draghi, C.; Avalos, J. B.; Rousseau, B. J. Chem. Phys. 2005, 122, 114503. (30) Packer, K. J.; Tomlinson, D. J. Trans. Faraday Soc. 1971, 67, 1302. (31) Be´e, M. Quasielastic Neutron Scattering; Hilger: Bristol, U.K., 1988. (32) Cabral, J. T.; Luzar, A.; Teixeira, J.; Bellissent-Funel, M. C. Physica B 2000, 276, 508. (33) Bordallo, H. N.; Herwig, K. W.; Luther, B. M.; Levinger, N. E. J. Chem. Phys. 2004, 121, 12457. (34) Teixeira, J.; Bellissent-Funel, M. C.; Chen, S. H.; Dianoux, A. J. Phys. ReV. A 1985, 31, 1913. (35) Di Cola, D.; Deriu, A.; Sampoli, M.; Torcini, A. J. Chem. Phys. 1996, 104, 4223. (36) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (37) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, U.K., 1987; Vol. 1987. (38) Verlet, L. Phys. ReV. 1967, 159, 98. (39) Evans, D. J.; Murad, S. Mol. Phys. 1977, 34, 327. (40) Sears, V. F. Can. J. Phys. 1967, 45, 237. (41) Egelstaff, P. A. An Introduction to the Liquid State; Academic: London, 1967. (42) Skaf, M. S. J. Chem. Phys. 1997, 107, 7996. (43) Skaf, M. S. Mol. Phys. 1997, 90, 25. (44) Svishchev, I. M.; Kusalik, P. G. J. Phys. Chem. 1994, 98, 728. (45) Laage, D.; Hynes, J. T. Science 2006, 311, 832. (46) Ropp, J.; Lawrence, C.; Farrar, T. C.; Skinner, J. L. J. Am. Chem. Soc. 2001, 123, 8047. (47) Svishchev, I. M.; Kusalik, P. G. J. Chem. Soc., Faraday Trans. 1994, 90, 1405. (48) van der Spoel, D.; van Maaren, P. J.; Berendsen, H. J. C. J. Chem. Phys. 1998, 108, 10220. (49) Lynden-Bell, R. M.; Steele, W. A. J. Phys. Chem. 1984, 88, 6514. (50) Ivanov, E. N. SoV. Phys. JETP 1964, 18, 1041. (51) Laage, D.; Hynes, J. T. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 11167.