Hydrogen Bonding and Dielectric Spectra of Ethylene Glycol-Water Mixtures from Molecular Dynamics Simulations.

Mixtures of ethylene glycol with water are a prominent example of media with variable viscosity. Classical molecular dynamics simulations at room temperature were performed for mixtures of ethylene glycol (EG) and water with EG mole fractions of xE = 0.0, 0.1, 0.2, 0.4, 0.6, 0.9, 1.0. The calculated dielectric loss spectra were in qualitative agreement with experiment. We found a slightly overestimated slowdown of the dynamics with increasing EG concentration compared to experimental data. Statistics of the hydrogen bond network and hydrogen bond lifetimes were derived from suitable time correlation functions and also show a slowdown in the dynamics with increasing xE. A similar picture is predicted for the time scales of EG conformer changes and for molecular reorientation. A slight blue shift was obtained for the power spectra of the molecular center of mass motion. The results were used to give a qualitative interpretation of the origin of three different relaxation times that appear in experimental complex dielectric spectra and of their change with xE.

Introduction

Water and ethylene glycol (EG) are perfectly miscible and the viscosity of EG is much larger than that of pure water. Therefore, EG–water mixtures are a prominent example of media with controllable viscosity. This makes them useful for investigations of solvent dynamics effects on electron transfer reactions (i.e., the dependence of the mechanism of elementary act on solvent viscosity, saddle point avoidance)[1,2] that are sometimes the dominant contributions to the reaction rate. The structural and dynamical details of water/EG mixtures at the molecular level are much more complicated than those of pure components. The interpretation often involves the properties of the hydrogen bond network in the mixtures. Experimentally, dielectric spectra of EG–water solutions have been reported in the classical region of frequencies (ℏω ≪ kB).[1] The imaginary part of the complex dielectric spectra was fitted to a Debye model with three characteristic relaxation times. The authors suggested a possible molecular origin of these modes at a very general and qualitative level in terms of H-bonding. Furthermore, H-bonding in EG–water mixtures was probed recently in the quantum (high frequency) region by Fourier transform infrared spectroscopy (FTIR),[3] broad-bandwidth sum frequency generation spectroscopy (BBSFG),[4] and near-infrared (NIR) technique.[5] According to the quantum mechanical theory of charge transfer, the region above 500–600 cm–1 can be considered the quantum (high frequency) region. From the viewpoint of modern quantum mechanical theories of electron transfer in condensed media,[6] classical solvent modes contribute to the activation (Franck–Condon) energy barrier, while their characteristic relaxation times are responsible in the adiabatic (strong coupling) limit for solvent dynamical effects. In calculations of corrections to kinetic parameters, quantum (i.e., high frequency) solvent modes usually reduce the activation barrier but might slightly decrease the rate constant of electron transfer due to the tunneling effect.[6]

A detailed understanding of the origins of the characteristic molecular relaxation times, in order to identify key modes in complex mixed solvents is of great fundamental interest but remains a complicated and challenging issue so far. In a previous article,[7] we provided a comprehensive analysis of H-bonding in liquid EG based on molecular dynamics (MD) simulations with a critical look on the differences between various force fields. On the basis of these findings we present new results for the more complicated case of EG–water mixtures. We also try to interpret the behavior of the imaginary part of the experimental dielectric spectra presented in ref (1).

This paper is organized as follows: Details of the MD simulations are reported in Section . Static and dynamic results of the simulated mixtures are discussed in Section , including dielectric spectra, hydrogen bond statistics and lifetimes, selected time correlation functions, and infrared power spectra. In Section we discuss possible relations of our results to the experimentally determined relaxation times.

Computational Details

The MD simulations were performed using the DL_POLY_4 package.[8] Mixtures of EG and H2O were prepared in cubic boxes. The mixtures were initially generated by randomly replacing water molecules with EG molecules, followed by an extensive equilibration. nE EG molecules and nW H2O molecules were mixed according to the mole fractions xE = nE/(nE + nW) = 0.0, 0.1, 0.2, 0.4, 0.6, 0.9, and 1.0. The respective numbers are given in Table .

Table 1

Mole Fraction of EG (xE), Number of EG (nE), Number of H2O (nW), Mass Fraction wE in Weight Percent, Average Box Length L in Å, Mean Density ρ̅ in kg/m3, Experimental Density ρexp,[9] and Static Dielectric Constants ε(0)a

xE	nE	nW	wE (%)	L	ρ̅	ρexp	ε(0)b	ε(0)c
0.0	0	503	0.00	24.70 ± 0.04	999.3	997.1	73.1	80.08
0.1	61	547	27.76	27.90 ± 0.03	1042.1	1035.3	77.4	72.71
0.2	122	486	46.38	29.35 ± 0.04	1071.0	1058.8	80.1	67.18
0.4	243	365	69.64	31.92 ± 0.04	1106.2	1084.9	80.0	57.6
0.6	365	243	83.81	34.13 ± 0.04	1126.7	1098.2	(87.7)	50.55
0.9	547	61	96.86	37.12 ± 0.04	1141.6	1108.6	(88.8)	43.69
1.0	512	0	100.00	35.87 ± 0.04	1146.6	1110.8	(79.1)	41.78

For pure SPC/E water a dielectric constant of 71.03 was reported in ref (10).

Values from this work; values in parentheses are poorly converged within 100 ns simulation time and are possibly subject to large errors.

Experimentally measured static dielectric permittivity.[1]

The mass fraction expressed in weight percent wE = nEmE/(nWmW + nEmE) is given in the fourth column of Table , where mE = 62.06784 u is the mass of EG and mW = 18.01528 u is the mass of H2O. After equilibration in the NVE ensemble followed by heating and cooling with the Nosé–Hoover thermostat,[11,12] final equilibration runs were carried out in the constant pressure ensemble at p = 1 atm and T = 298 K with a thermostat relaxation time of 1 ps and a barostat relaxation time of 3 ps, using a Berendsen barostat and thermostat (NpT).[13] The same conditions were applied to the subsequent production runs. A more detailed description of the equilibration procedure is included in our previous article dealing with pure liquid EG.[7] We are aware of the deficiencies of the Berendsen ensemble with respect to the sampling of canonical trajectories. Nevertheless, the Berendsen ensemble is widely used and comes close to experimental conditions. In our ensemble the volume of the simulation box and the density change slightly in each time step to keep the pressure constant. The mean values of the densities are also given in Table . They are close to the interpolated experimental values ρexp.[9] For example, pure EG has an experimental density of 1.11 g/cm3 at 20 °C,[14] while our value is 1.15 g/cm3 at 25 °C.

In our previous work on pure EG we carefully tested different force fields of the OPLS type and generally obtained good agreement for radial distribution functions but several rather unsettling results for the crucial dihedral angle distributions. The OPLS-AA-SEI-M force field[15] with scaled electrostatic interactions (scaling 1–4 van der Waals interactions with 0.65 and 1–4 Coulomb with 0.6, 1–5 interactions with 0.8) yielded dihedral distributions in acceptable agreement with experimental data. It is a successor of the OPLS-AA-SEI[16] parameter set that in turn is based on the original all atom OPLS-AA force field.[17] The importance of scaling factors was discussed by Lin et al., who found that a 1–4 scaling of the Coulomb forces between 0.6 and 0.8 gives best results for diffusion coefficients, the thermal conductivity, and the H-bond network.[18]

The OPLS-AA-SEI-M force field for EG was used together with SPC/E for water[19] for all simulations in the present work. The SPC/E water model was chosen because it is especially accurate for capturing experimental properties of water such as the diffusion coefficient and dielectric constant, and our work relates strongly to these properties. It also performed well in combination with OPLS-AA for hydration of amino acid analogues.[20] The combination rules for OPLS force fields were used to define the parameters between EG and H2O, . Periodic boundary conditions were employed. We used Ewald summation to handle the long-range electrostatic forces and a cutoff of 10 Å for short-range forces. Depending on the quantity of interest, we selected an integration time step of 1 or 2 fs with simulation times T from 100 ps up to 100 ns. The recording strides varied in order to make efficient use of computing resources. The Velocity Verlet scheme was used to integrate the equations of motion.

Results and Discussion

Dielectric Loss Spectra and Molecular Diffusion

The complex dielectric function ϵ(ω) can be estimated without applying external electric fields by using linear response theory using the formula derived by Neumann et al.[21] for MD simulations with periodic boundary conditions and Ewald summation. Together with conducting boundary conditions ϵRF = ∞ outside the Ewald repetition units the dielectric spectrum is given bywhereis the autocorrelation function of the total dipole moment M⃗(t) of the whole simulation box, Text = 298 K here, L is the Laplace transform, ⟨M2⟩ is the mean squared dipole moment, and V is the box volume. The total dipole moment is calculated as a sum over dipole moments of individual molecules

Molecules that are split at the boundaries due to the periodic boundary conditions have to be reconstructed before calculating p⃗(t). The static dielectric constant is given by the equation

In the NpT ensemble, V fluctuates slightly with time. In eq we used its mean value V = ⟨V(t)⟩. We also performed a direct comparison between the Berendsen NpT ensemble and the Nosé–Hoover NVT ensemble and found significant differences in neither ϕ(t) nor ϵ(ω).

With increasing EG content the dynamics slow down drastically and much longer simulation times are needed for convergence of the dipole correlation functions. The longest simulations we could perform reached up to 100 ns, which is still not long enough to reach complete convergence for EG mole fractions equal to or larger than xE = 0.6. For these long simulation times we modified the DL_POLY code by calculating M⃗(t) at each time step on the master node. The time-consuming writing to the history file was suppressed. This achieved a speedup factor of 5–10, still requiring a total simulation time of about one month on a dedicated single 8 core Xeon X5650 CPU.

The imaginary parts of the dielectric spectra calculated for selected values of the EG content are displayed in Figure . Note that the tail of the dipole correlation function ϕ(t) was replaced by a fitted single exponential function; this procedure smooths the curve but has only a negligible effect on the positions and heights of the maxima. We refer to ref (22) for a discussion about the spectra for pure H2O with the SPC/E force field and their comparison with experiments. The xE dependencies of ϵ look qualitatively similar to those obtained experimentally in ref (1), since all calculated spectra exhibit a strong maximum which shifts smoothly to lower frequencies with increasing EG concentration. However, the shift of the maxima of the computed spectra to the red region is overestimated compared to the experimental data and the decreasing value of the absorption maximum is not reproduced by our model. This result and the deviations for xE > 0 in ϵ(0) reported in Table 1 in the Supporting Information may be caused by the fact that the molecular dipole moment of EG differs between our model and accurate quantum chemical data. However, it should be noted that only a weak correlation between the molecular dipole and ϵ(0) has been found in the case of several water force fields.[23] The average value of 3.9 D in our simulation of pure EG is larger than the literature value of 2.34 D[24,25] for the liquid and also larger than the conformationally weighted dipole moment of 2.26 D[26] from quantum chemical calculations at the CCSD(T)/aug-cc-pVTZ level. A scaling of the atomic charges of the EG model toward these presumably more accurate values would influence all other carefully adjusted properties of the OPLS-AA-SEI-M force field and was therefore not attempted. Nevertheless, the qualitatively correct results provide a plausible picture on the interplay of different solvent modes.

Figure 1

Calculated dielectric loss spectra for EG–water mixtures for different EG mole fractions xE.

The diffusion coefficients of EG and water molecules were calculated both via the velocity autocorrelation function and from the Einstein equation. The velocity autocorrelation functions of the center of mass velocitiesare related to the self-diffusion constants according to

This method was applied to the 100 ps simulations with a recording stride of 5 fs. The Einstein relationwas employed for simulations with total lengths of 1 ns and recording strides of 50 fs. Overall, for our systems the Einstein approach yielded better convergence and therefore more reliable results. The reason is that in the approach via the velocity autocorrelation function, D depends crucially on the long-time part of the velocity autocorrelation function (eq ), which in our case, especially for large values of xE, is quite noisy. We investigated possibilities to reduce the noise without reducing the information content of the data. It turned out that different approaches to calculate the correlation function do not influence it and it does not result from a coupling between thermo- and barostat, which is sometimes problematic. We have to conclude that it is intrinsic to our data and may only be removed by greatly increasing the simulation size and time. The results from both approaches are collected in Table . The self-diffusion coefficient of (2.8 × 10–9) ± (0.06 × 10–9) m2/s for SPC/E water, as given in ref (27) at a temperature of 298.2 ± 1.4 K, is slightly larger than our values of 2.5 × 10–9 m2/s (eq ) or 2.4 × 10–9 m2/s (eq ) at 298 K. An experimental value of DH2O = 2.299 × 10–9 m2/s was obtained at T = 25 °C from NMR measurements.[28] The self-diffusivity DEG = 9.61 × 10–11 m2/s, as measured for pure ethylene glycol, is 1 order of magnitude larger than our calculated value DEG = 1.2 × 10–11 m2/s and, again, as in the case of the dielectric spectra, points to an overstabilization of the hydrogen bond network by the force field. Rodnikova et al. reported an experimental value of DEG = 0.834 × 10–10 m2/s at 298 K from a spin echo experiment.[29]

Table 2

Calculated Self-Diffusion Coefficients DEG and DH2O for Several EG Concentrations xE in Mixed Solutions

	DH2O (m2/s)		DEG (m2/s)
xE	eq 7	eq 6	eq 7	eq 6
0.0 (pure H2O)	2.4 × 10–9 (3.8 × 10–11)a	2.5 × 10–9	-	-
0.1	1.4 × 10–9 (3.5 × 10–11)	1.5 × 10–9	5.2 × 10–10 (3.9 × 10–11)	5.5 × 10–10
0.2	8.6 × 10–10 (2.9 × 10–11)	8.6 × 10–10	3.1 × 10–10 (4.7 × 10–11)	3.6 × 10–10
0.4	3.5 × 10–10 (1.9 × 10–12)	3.6 × 10–10	1.2 × 10–10 (1.2 × 10–11)	1.2 × 10–10
0.6	1.2 × 10–10 (1.7 × 10–12)	1.7 × 10–10	3.5 × 10–11 (1.8 × 10–12)	5.7 × 10–11
0.9	4.3 × 10–11 (1.5 × 10–11)	9.8 × 10–11	1.6 × 10–11 (5.4 × 10–13)	3.1 × 10–11
1.0 (pure EG)	-	-	1.2 × 10–11 (9.6 × 10–13)	2.4 × 10–11

Estimated numerical errors are given in parentheses from the standard deviation of the self-diffusivities in x, y, and z directions. Possible systematic errors in the force field are not included by this procedure.

As can be seen from Table , the molecular mobility of both EG and water decreases significantly with increasing EG content in the mixture. We see in Figure 1S in the Supporting Information a nonlinear dependence of the self-diffusion coefficients on xE, which can also be seen in ref (30) for simulations with GROMOS-compatible force fields with and without inclusion of polarization (by a charge on spring model) for EG; both DEG and DH2O drop most significantly in the interval from xE = 0 to xE = 0.4.

It would be interesting if this feature could be explained in terms of the Einstein equationwhere m is the mass of a particle and γ is the medium friction coefficient.[31] It is known that the so-called “slow” and “medium” solvent relaxation times (τslow and τmid, respectively) contribute mainly to the friction coefficient of the solvent (γ = τω2, where ω is a characteristic solvent frequency).[31] Both τslow and τmid, however, increase monotonously and practically linearly in the region from xE = 0.0 to xE = 0.4.[1] The nonlinear behavior of D vs EG content could be described via an effective particle mass that depends on xE (eq ), which would contain information about the formation of strongly bound molecular associates with H-bond networks formed by the water molecules.

Time Scales of Molecular Motion from Correlation Functions

We used the reorientation correlation function to quantify the time scale of molecular reorientation. Given a unit vector u⃗(t) it is defined bywhere, for the sake of comparisons with NMR measure-ments, we use the second-order Legendre Polynomial P2(x) = (3x2 – 1)/2 to define the orientational metric. The average ⟨...⟩ goes over all times τ in the interval [0, T – t]. We also define a half-life time t1/2 by the equation to quantify the reorientation time. This definition of the half-life time is chosen throughout the manuscript instead of, for example, fitting exponential decay rates, in order to be able to compare the results of different correlation functions directly with each other. It is also preferred to define τ by integration of time-correlation functions because these are often poorly converged for large t. The vector u⃗ was chosen to be a unit vector along the C–C backbone axis for EG, and along the C2 symmetry axis for H2O. The results are shown in Figure . For both EG and H2O the reorientation slows down with increasing EG content even though H2O molecules naturally reorient much faster than EG molecules.

Figure 2

Dependence of the orientational lifetimes on the mixing ratio, as calculated from the reorientation correlation function of the EG backbone (dashed line connecting C–C) and the H2O symmetry axis (dashed line through O).

EG molecules can be classified roughly into gauche and trans states depending on the dihedral angles of the OCCO backbone and the two HOCC angles. For example, a gTt molecule would have the first HOCC angle gauche (g), the OCCO backbone trans (T), and the second HOCC angle also trans (t).[32] Trans–gauche transitions occur in our simulations at room temperature. The dynamics of conformer changes, i.e., the lifetime of the conformer flips, is of particular interest. To measure the state of a dihedral angle α at a time t we introduce a function fgt(α,t) defined aswhere αc is a critical dihedral angle that divides the angular distribution into gauche and trans subsets, where α < αc corresponds to the gauche state. Inspection of the dihedral angle distributions of pure EG[7] leads us to choose αc = 130° for the HOCC angle and αc = 120° for the OCCO angle. In the correlation functionthe averaging goes over all EG molecules and initial times τ in [0,T–t]. Finally, a normalization factor ensures xconf(0) = 1. The decay time of xconf is a measure for the time scale of conformer changes. To stay consistent, we again use the half-lifetime t1/2 defined as xconf(t1/2) = 1/2 as a measure of the flipping time. Figure shows the lifetimes from 1 ns simulation runs. The conformer changes slow down with EG content in a fashion similar to what has been observed above for the molecular reorientation. Among them, OCCO backbone conformer transitions occur with a higher frequency than the HOCC transitions. Figure , Figure , and Figure indicate a slowdown of the dynamics in the liquid with higher EG content, accompanied by a stabilization of the H-bond network. This network stabilization at room temperature is somewhat counterintuitive to the xE-dependence of the melting point.[33] The experimental melting point of EG–water mixtures drops from 0 °C for pure H2O toward −49.9 °C for 58 wt% (xE = 0.29). It slowly increases again via a small maximum of −42.6 °C at 77.5 wt% (xE = 0.50) to −43.4 °C for 81.1 wt% (xE = 0.55) and then steeply climbs to −12 °C for pure EG.[33] The intriguing region between the two eutectic points xE = 0.29 and xE = 0.55 contains our investigated mole fraction xE = 0.4 that also stands out in many properties investigated in the present work. For example, the excess molar volume of an EG–H2O mixture also displays a maximum at ca. xE = 0.4.[3] In both Figures and 3, we obtain a change in slope at xE = 0.4.

Figure 3

Lifetimes of the OCCO and HOCC dihedral angles defining the EG conformers.

Hydrogen Bonding

We identified hydrogen bonds using a geometrical criterion by using the procedure that we implemented before for pure ethylene glycol.[7] A donor is H-bonded to an oxygen acceptor if the O–H (ROH) and O–O (ROO) distances of the two participating hydroxyl groups are smaller than critical distances given by the first minima of the corresponding radial distribution functions (rdf) and the H–O···O valence angle ϕ is smaller than 30° (intermolecular H-bonds) or 60° (intramolecular H-bonds). Complexity increases compared to pure EG because four types of H-bonds exist in mixtures EGA – H2OD, H2OA – EGD, EGA – EGD, H2OA – H2OD, with A for acceptor and D for donor. Critical ROO distances for EGA – EGD are in the range of 3.1–3.26 Å, and for H2OA – H2OD as well as mixed bonds, ROO = 3.26 to 3.4 Å. Critical ROH distances are found between 2.3 and 2.4 Å with EGA – EGD and EGA – H2OD, varying strongly with concentration, whereas for H2OA – H2OD and H2OA – EGD they remain rather constant around 2.33 Å. The ROH and ROO dependences on xE (not shown here) indicate that EGA – EGD H-bonds are strongest at xE = 0.4 (smallest ROH), and that H2OA – H2OD bonds are the weakest ones (largest ROO).

For all types of H-bonds, the calculated angular distributions of cos(ϕ) clearly favor linear H-bonds as compared to a uniform distribution of the two spatial vectors. These distributions are exemplified for xE = 0.4 in Figure . However, in absolute numbers, the majority of the H-bonds has angles between 5° and 15° and angles much smaller than 5° are unlikely because of the small available solid angle.

Figure 4

Distribution of the cosines of O–H–O angles in different H-bonds at xE = 0.4.

On average, a molecule is involved in nHB intermolecular H-bonds; The values of nHB are shown in Figure for EG and for H2O molecules, where the contributions from different types of H-bonds are color-coded. For pure H2O nHB(H2O) = 3.42 and for pure EG nHB(EG) = 3.86; Padró et al. predicted nHB(EG) = 3.9.[34] Jorgensen et al. obtained values of nHB (H2O) = 3.5 to nHB (H2O) = 3.73 in simulations of pure water for different force fields, not including SPC/E.[35] In addition to the contributions to nHB(EG) and nHB(H2O) from the four types of H-bonds, we also show the average number of H-bonds, nHB(EG, H2O), per either of the two molecules in both panels in Figure (black circles). nHB(EG, H2O) × (nE + nW)/2 is equal to the average of the total number of H-bonds in the simulation box. nHB(EG, H2O) increases continuously from pure water to pure EG. For H2O molecules mixed bonds (EGA – H2OD, H2OA – EGD) dominate over like bonds (H2OA – H2OD) for xE > 0.4. For EG molecules, mixed bonds dominate over like bonds (EGA – EGD) for xE ≤ 0.4. Note that the xE resolution is limited by the small number of data points. EG is the preferred acceptor molecule for mixed bonds. For low EG concentrations, unexpectedly large values for nHB(EG) were obtained with 4.79 H-bonds per EG molecule in the xE = 0.1 mixture; in most cases (51.4%) water donates hydrogen to oxygen of EG.

Figure 5

Average number of H-bonds per H2O molecule, nHB(H2O) (left), and per EG molecule, nHB(EG) (right), and their decomposition into different types. The black bullets give the average number of H-bonds per either one of the two molecules, nHB(EG, H2O), (XA–XD).

Intramolecular H-bonds are rare, at least in our simulations with the OPLS-AA-SEI-M force field and the present H-bond criterion.[7] They seem to occur more frequently in simulations with the AMBER99 force field where their number presumably increases with temperature.[36] Their real amount is unclear at present. No intramolecular hydrogen bonds were found in a near-infrared spectroscopy study of EG–H2O mixtures.[5] At high pressure above 4 GPa a liquid–solid transition was observed by changes in the hydrogen bond network caused presumably by formation of intramolecular H-bonds.[37] In a recent spectroscopy study by Kollipost et al.[38] a highly symmetric dimer structure with four intermolecular and no intramolecular H-bonds was dominant and most stable, which points to a minor role of intramolecular H-bonding in contrast to the monomer. This gives rise to the assumption that the molecular geometry in pure and mixed liquids can also be very different from the most stable monomer geometries found in quantum chemical calculations. Interestingly, dimer structures stabilized by intramolecular H-bonds have been proposed in solution in an earlier study.[39] However, both studies agree that EG prefers the gauche conformation in dimers in both gas phase and liquids. Other MD simulations of pure EG also indicate that the concentration of dimers in the liquid is rather low.[36]

Individual lifetimes of H-bonds and population correlation functions SHB(t, t*) were evaluated for the four types of H-bonds and concentrations.[7,40] We used again the half-life time t1/2 defined by SHB(t1/2, t*) = 1/2 without exponential fitting as a quantitative measure. An overview of t1/2 is shown in Figure for t* = 0.05 ps in the left panel and for t* = 5 ps in the right panel. We observe that the lifetimes of all H-bond types increase with EG concentration. This is not surprising for EGA – EGD bonds since EG is a heavier molecule than H2O and also indicates that EGA – EGD bonds are stronger than H2OA – H2OD bonds (see, for example, ref (41)) which is also reflected in a slightly smaller average H-bond length. It is also the case, though, for EG–H2O and H2O–H2O H-bonds, indicating a stabilization of the whole H-bond network. The lifetimes are generally much smaller for small t* = 0.05 ps than for large t* = 5 ps. This indicates that most of the bonds break regularly over short time intervals and are formed again, whereas the time scales of irreversible destruction by diffusion are longer. EGA – EGD bonds live longer than H2OA – H2OD bonds and mixed bonds live longer if EG is the donor. Overall, the EG H-bond network is stronger and more long-lived than the H2O network. At xE = 0.4, we see a slight bend in the EGA – EGD and H2OA – EGD lifetimes for large t*. Belashchenko et al. calculated a hydrogen-bond lifetime of 147 ps for pure EG at 323 K with a different method, where only the complete diffusive H-bond destruction was considered.[36]

Figure 6

Hydrogen-bond lifetimes defined by SHB(t1/2, t*) = 1/2, for strides t* = 0.05 ps (left) and t* = 5 ps (right).

After completion of this study, work has been published[42] in which EG/water mixtures have been simulated both at atmospheric pressure and at 100 MPa. The hydrogen bonding was studied by an approach similar to the one in this work. There, an EG force field was used[43] which we had investigated previously but did not employ due to its deficiencies in the dihedral angle distributions.[7] The dynamic results concerning H-bonding (lifetimes) in this study are similar to our findings. It is, however, surprising that in ref (42) the average number of intermolecular hydrogen bonds per molecule decreases drastically from pure water to pure EG, in contrast to what we showed in Figure , and water–water H-bonds are predicted to have a longer lifetime than EG–EG H-bonds, although the overall dynamics of the H-bond network slow down with increasing xE, similar to our results given in Figure .

Power Spectra

Power spectra were calculated from normalized velocity autocorrelation functions (Equation ) by taking the real part of their Fourier transforms. The power spectrum of the center of mass motion is shown for EG molecules in Figure and for H2O molecules in Figure . Experimentally determined Rayleigh/Raman spectra also show a maximum, however, slightly above 50 cm–1 for pure EG and with a different line shape.[44] A peak at 50 cm–1 corresponds to a time scale of 0.7 ps.

Figure 7

Power spectra of the EG molecules, obtained from the velocity autocorrelation functions of the centers of mass of the EG molecules.

Figure 8

Power spectra of H2O molecules, obtained from the velocity autocorrelation functions of the centers of mass of the H2O molecules.

The wide shoulder in the power spectra of water molecules was attributed by Balucani et al. to both longitudinal and transversal contributions.[45] This shoulder becomes weaker for larger xE. The power spectra shown in Figure (EG) and Figure (water) reveal a slight blue shift of the maxima with increasing EG fraction. While all previously discussed dynamical quantities slow down with increasing EG content, the intra/intermolecular femtosecond dynamics get slightly faster. The decrease of the self-diffusion coefficient for larger xE, already reported in Table , is also seen from the decrease of intensities at ν = 0 cm–1 in Figures and 8.

Interpretation of Solvent Relaxation Times

In this section we discuss a possible origin of the solvent relaxation times that appear in the dielectric response of EG–water mixtures. These three times are reported in ref (1) where they were derived from experimental data and revealed qualitatively interesting behavior. The first (slowest, <150 ps) time exhibits a practically linear dependence with xE. The plot of the second (medium, <35 ps) relaxation time versus xE has a plateau-like maximum and an inflection point at xE = 0.4. The third (fast) time ranges from 2 to 4 ps; however, these values have large error bars and the shape of τfast(xE), presented in ref (1), seems hardly reliable. Most of the experimental dielectric spectra reported in the literature describe basically classical (low frequency) solvent modes. Quantum (high frequency) modes can be addressed as well, but this requires additional measurements using vibrational spectroscopy (see, for example, ref (46)).

Since in our computational study the lifetimes obtained from the EG reorientation correlation function were found to be the slowest (Figure ), this quantity might be responsible for the origin of the slowest relaxation time in the experimental dielectric response. In turn, it can be seen that the scale of the calculated lifetimes of H2O–EG and H2O–H2O H-bonds versus xE, as shown in Figure , is comparable with that of the medium relaxation time.[1] Therefore, we can argue that this relaxation time stems from the dynamics of different H-bonds in water–EG mixtures.

The third (fast) relaxation time of dielectric response might in turn be attributed to some low-frequency vibration modes in EG–H2O and EG–EG associates (Figure , Figure ). To reinforce this statement, we calculated the harmonic vibrational frequency spectra of selected EG–water and EG–EG associates by DFT (see details in the Supporting Information). Some typical snapshots of hydrogen-bonded molecular clusters with large lifetimes were initially extracted from the MD-data. Then, the geometries of these snapshots were used as the starting geometry of associates in DFT calculations. Characteristic times τeff of low-frequency intermolecular vibration modes in the far-infrared region can be defined as reciprocal frequency values. As seen from Figure 2S, the estimated τeff range from 0.5 to 2.6 ps and indeed fall into the region of the experimentally detected fast relaxation times.[1] Thus, our results also shed light on the molecular origin of dynamic properties of mixed solvents, at least qualitatively. A similar approach could be extended to other solvent mixtures which are used as media with variable viscosity (for example, water–glycerol).

Conclusions

We simulated mixtures of ethylene glycol and water with concentrations varying from pure H2O to pure EG using the OPLS-AA-SEI-M force field together with SPC/E water. Dielectric spectra were calculated from very long simulations and show a qualitative agreement with experiments with respect to the red-shift of the pronounced maximum in the GHz regime with increasing EG content. This maximum shifts from 15 GHz for pure H2O to 6.2 GHz for the mixture with an EG mole fraction xE = 0.1 and further to 0.4 GHz for xE = 0.6, whereas in the experiment it shifts from 9 GHz for xE = 0.1 to 2 GHz for xE = 0.6. We also obtained a marked slowdown of the dynamics of the hydrogen bond network, of the molecular reorientation, and of the flipping times between EG conformers. This slowdown was accompanied by an overestimated (compared to experiments) decrease of the self-diffusion coefficients. A slight blue shift was found in the femtosecond dynamics from calculated power spectra of the molecular center of mass motion. An attempt was made to explain the origin of three different relaxation times in the experimentally observed dielectric spectra. The slowest relaxation time (<150 ps) is presumably attributed to the reorientation of EG molecules, the medium one (<35 ps) might be related to diffusion-limited (large t*) hydrogen-bond lifetimes, and the fast one (<4 ps) could result either from fast breaking and rebonding of hydrogen-bonds due to (collective) librations/vibrations (t* small) or from low-frequency vibrations of EG clusters or molecules. The statistics of the hydrogen bond network showed a monotonous increase in the average number of H-bonds per molecule from 3.42 for pure H2O to 3.86 for pure EG. Up to 4.79 H-bonds per EG molecule were found in the dilute xE = 0.1 mixture. Although some of our findings are obviously of qualitative nature, the results provide a useful insight into the complex dynamics of EG–water mixtures. Further improvements in force fields for simulations of EG–water binary liquids seem necessary for quantitative predictions of the dynamics. In particular, the average dipole moment of EG molecules in the liquid could be improved.