Theoretical Study of the Temperature Dependence of the Vibrational

Jul 22, 2014 - variation of the hydrogen bonds network in the liquid. The vibrational bending ... molecule, to the degrees of freedom of the network f...
0 downloads 0 Views 702KB Size
Article pubs.acs.org/JPCB

Theoretical Study of the Temperature Dependence of the Vibrational Relaxation of the H2O Bend Fundamental in Liquid Water and the Subsequent Distortion of the Hydrogen Bond Network Beatriz Miguel,† José Zúñiga,‡ Alberto Requena,‡ and Adolfo Bastida*,‡ †

Departamento de Ingeniería Química y Ambiental, Universidad Politécnica de Cartagena, 30203 Cartagena, Spain Departamento de Química Física, Universidad de Murcia, 30100 Murcia, Spain



ABSTRACT: The molecular dynamics with quantum transitions method is used to study the temperature dependence of the relaxation dynamics of the H2O bend fundamental in liquid water in the range from 277 to 348 K and the subsequent variation of the hydrogen bonds network in the liquid. The vibrational bending degrees of freedom of the water molecules are all described by quantum mechanics while the remaining translational and rotational motions are described classically. The participation of the H-bonds in the relaxation process is studied taking into account the dependence of the relaxation lifetimes on the number of H-bonds formed by the initially excited water molecule and the amount of energy transferred into the hindered rotational motions. It is found that the intermolecular vibrational energy transfer plays an important role in the relaxation mechanism, with almost no temperature dependence, and that the energy transfer into the rotational degrees of freedom is favored over the energy transfer into the translational motions. The thermalization of the system after relaxation is completed in a time scale shorter than the time taken for the H-bond network to recover. The relaxation and equilibration times calculated compare well with experimental and previous theoretical results.

1. INTRODUCTION The use of the temperature as a variable in ultrafast spectroscopic techniques has provided new insights into the study of the vibrational relaxation processes that occur in bulk water.1−11 The first experimental studies that considered the temperature dependence of the vibrational lifetimes in water were carried out by Bakker et al.1,2 and dealt with the relaxation of the OH stretch of the HOD molecule dissolved in liquid D2O. It was observed that the OH stretch relaxation lifetimes increased from 0.74 ± 0.01 to 0.90 ± 0.02 ps when raising the temperature from 298 to 363 K. This lifetime increase was somewhat unexpected since simple perturbative arguments indicated that the coupling between the initially excited vibrational state and the remaining states of the system should become stronger with thermal agitation and thus cause an increase of the relaxation rate. Bakker et al. suggested that this anomalous temperature dependence was due to the role played by the hydrogen bonds (H-bonds) as dominant accepting modes of the excess energy released by the OH stretch. The rise of the temperature would accordingly weaken the H-bonds and consequently results in an increase of the relaxation time. Bakker et al. later conducted similar experiments in the pure liquid water3 (H2O) showing that the OH stretch also relaxes slower as temperature rises, with measured lifetime values increasing from 260 ± 18 fs at 258 K to 320 ± 18 fs at 358 K. In this case, however, the overtone of the bending mode (2δH2O) was argued to be the main energy acceptor during the relaxation process and the responsable of the anomalous © 2014 American Chemical Society

temperature dependence. On raising the temperature, the frequencies of the OH stretch and the bending modes would be respectively blue- and red-shifted, resulting in an increase of the energy gap that would eventually slow down the relaxation process. Although this mechanism should also work in the HOD/D2O system,12,13 further experiments4 showed that the excited OH stretch state decays back in this case to the ground vibrational level without appreciably populating any intermediate states. The relaxation lifetime of the OD stretch of HOD in liquid H2O also exhibits the anomalous increase with temperature,7−10 and a direct v = 1 → 0 relaxation mechanism of the excited OD stretch has also been proposed in this case to account for it.8 In addition, Tokmakoff et al.9,10 have used temperature-dependent two-dimensional infrared spectroscopy of this system to investigate H-bond arrangements, showing that the time scales of spectral diffusion and reorientational relaxation decrease from 2.4 to 0.7 ps and 4.6 to 1.2 ps, respectively, over the temperature range between 278 and 345 K. The difficulties found in unraveling the relaxation pathways of the high-frequency stretching modes have favored the study of the relaxation of the water bending mode.11,14−20 The lower frequency of this mode dramatically simplifies the relaxation mechanism, since the excess vibrational energy of the bend Received: June 12, 2014 Revised: July 22, 2014 Published: July 22, 2014 9427

dx.doi.org/10.1021/jp5058447 | J. Phys. Chem. B 2014, 118, 9427−9437

The Journal of Physical Chemistry B

Article

system, and the computational details of the simulations. The numerical results are presented, discussed, and compared with previous theoretical treatments and with experiments in section 3, and conclusions are given in section 4.

fundamental can only be transferred to the rotation of the water molecule, to the degrees of freedom of the network formed by the surrounding water molecules or to the bending mode of another water molecule through intermolecular vibrational-tovibrational (VV) energy transfer, as observed for the OH stretching modes.21,22 Elsaesser et al.14,15 experimentally studied the relaxation of bend fundamental of pure H2O by femtosecond pump−probe spectroscopy finding a lifetime value of 170 fs at room temperature, which was later confirmed in subsequent experiments.16,17 This short lifetime is determined by the coupling of the solute molecule with the fluctuating liquid environment forming the H-bond network. It was also found that the energy transfer to the solvent induced a red-shift and a slight narrowing of the bending absorption characterized by a time constant of 0.77 ps, which accounts for the thermalization time scale of the system. The anisotropy decay measured in the Elsaesser et al.14 experiments was identical to the depopulation of the bend fundamental, showing a rapid delocalization of vibrational energy over a large number of water molecules. More recently, Ashihara et al.11 have carried out mid-infrared pump−probe spectroscopy measurements of the bending mode relaxation in liquid H2O at different temperatures, showing that the lifetime increases from 170 ± 15 fs at 295 K to 250 ± 15 fs at 348 K. This behavior was explained by the decrease in the spectral overlap between the bending and the high-frequency librational modes, which constitute the dominant relaxation pathway. Interestingly, these authors argued that previous polarization-resolved anisotropy experiments15 indicated a constant pump−probe anisotropy during the bend lifetime, ruling out any important contribution from the intermolecular bend-to-bend vibrational energy transfer route. From a theoretical point of view, the simulation of the dependence of the vibrational relaxation of liquid water with temperature turns out to be more challenging, as shown by the comparatively fewer number of studies of this type performed.13,23−25 In a first pioneer work, Lawrence and Skinner13 analyzed the temperature dependence of the OH stretch lifetime of HOD in liquid D2O using a hybrid quantum/ classical perturbative approach. More recently, Hynes et al.23 studied the temperature dependence of the water reorientation using molecular dynamics (MD) simulations and later extended25 their classical methodology to track preliminarily the energy transfer in the vibrational bending relaxation in water, and Paesani24 investigated the temperature dependence of the linear and nonlinear infrared spectra of HOD dissolved in H2O through the application of a first-principles approach. In this work, we have investigated the temperature dependence of the relaxation dynamics of the bending mode in liquid water and the subsequent distortion of the hydrogen bond network in the liquid using the molecular dynamics with quantum transitions (MDQT) method. This method was originally developed by Tully26,27 to describe electronic transitions and has been adapted by our group to simulate vibrational transitions.28−32 As in our previous work on the relaxation of the HOD bend in liquid D2O,32 we describe all of the bending modes of the H 2 O molecules quantum mechanically in order to treat the VV transitions naturally and properly quantify the relative importance of the intermolecular energy transfer processes as well as the role played by the H-bond network in the relaxation pathways. The paper is organized as follows. In section 2 we describe the vibrational MDQT method, its application to the liquid water

2. METHODOLOGY 2.1. The Vibrational MDQT Method. The vibrational MDQT method is an adaptation of the original electronic MDQT method proposed by Tully26,27 to simulate vibrational transitions, and it has proved to be very useful in the description of vibrational relaxation processes.29,32−34 The method has been developed and presented extensively in our previous works,32 so we just outline its main features here. As in any hybrid quantum/classical method,35−42 the system under study is divided into a classical subsystem described using time-dependent coordinates and conjugate momenta and a quantum subsystem characterized by a time-dependent wave function depending on the quantum coordinates. The Hamiltonian of the entire system is then expressed as a sum of the Hamiltonians of the quantum and classical subsystems plus a term containing the kinetic and potential couplings between both subsystems. The forces acting on the classical subsystem are calculated by averaging the total Hamiltonian over the wave function of an individual state of the quantum subsystem, and simultaneously, the time-dependent wave function that describes the quantum subsystem is propagated using the Hamiltonian of the system evaluated at the present values of the classical coordinates. The quantum equation of motion is usually solved by expanding the time-dependent wave function in a given basis set. Although the best option in principle is to use an adiabatic basis set, the need to evaluate it at each time interval supposes a high computational cost which may become prohibitive in systems like the present one in which a large number of quantum states participate in the vibrational relaxation. The alternative is then to use diabatic basis sets, which do not depend on the classical coordinates and therefore remain unchanged during the quantum propagation.43 To reinforce this approach, we should note that in previous studies28−34 on the vibrational predissociation of van der Waals clusters and the vibrational relaxation of liquid water the use of diabatic vibrational functions has been shown to provide results for different magnitudes that compare well with those measured experimentally and also with values calculated using different theoretical approximations, in contrast with the poor performance usually offered by the diabatic basis in electronic MDQT calculations. In the MDQT method,26 the transition probabilities between the quantum state that governs the motion of the classical subsystem and the remaining quantum states are calculated in such a way that the fraction of classical trajectories in the ith quantum state, the so-called classical population, is equal at every time to the quantum populations. Also in the method, when a hop between two quantum states is invoked, the classical momenta are adjusted in order to keep the total energy of the system constant.26,44,45 As in our previous work32 about the vibrational relaxation of the HOD bend fundamental in liquid D2O, the Hamiltonian for the liquid H2O system is written in terms of normal coordinates to describe the vibrational motion of each molecule, Euler angles for the rotational motions, and the center of mass vectors for the translational motions. Since the frequencies of the symmetrical and antisymmetrical stretching modes of the H2O molecules are much higher than the frequency of the 9428

dx.doi.org/10.1021/jp5058447 | J. Phys. Chem. B 2014, 118, 9427−9437

The Journal of Physical Chemistry B

Article

Table 1. Average Number of H-Bonds (nHB), Number of H2O Molecules in the First (nfss) and Second (nsss) Solvation Shells, Percentage and Lifetimes of the H-Bond Environments and H-Bond Lifetimes at Different Temperatures Obtained in Hybrid Quantum/Classical Equilibrium Simulations T (K) ρ (g/cm3)a nHB nfss nsss H-bond environment (%)

H-bond environment lifetime (fs)

H-bond lifetime (fs) a

1−2 3 4 5−6 1−2 3 4 5−6 mean value

277 0.999975 3.83 4.29 18.51 2.94 20.8 66.4 9.86 43 62 114 40 78 557

295 0.997773 3.78 4.32 18.21 4.33 24.3 60.5 10.9 41 59 94 36 67 488

323 0.98804 3.70 4.30 17.98 6.52 28.2 54.0 11.2 42 52 72 31 55 386

348 0.97484 3.63 4.26 17.66 8.77 31.1 49.0 11.1 37 47 60 27 47 328

Reference 60.

bending mode,46 only the bend fundamentals are expected to participate in the relaxation process, and accordingly, we freeze the stretching normal modes at the molecular equilibrium geometries and approximate the bending normal modes as harmonic oscillators. All the simulations performed include the rovibrational couplings in the Hamiltonian of every water molecule, which nevertheless have been shown32 to play a minor role in the energy redistribution that follows the vibrational relaxation of the bend fundamental in the HOD/ D2O system. Details of the use of both the diabatic basis set employed to expand the time-dependent wave function of the quantum subsystem and the second-order expansion of the intermolecular potential around the molecular equilibrium geometries to evaluate the potential energy matrix elements have also been given in our previous study on the HOD/D2O system.32 We have conducted two different sets of simulations. In the first, which we refer to as MDQT-A, the quantum subsystem is described using Nw + 1 vibrational states (with Nw being the number of water molecules in the simulation), each expressed as a product of Nw one-dimensional harmonic bending wave functions. The ground vibrational state is composed of the v = 0 states of every molecule, while the remaining Nw excited vibrational states are built upon one quantum bending excitation of a different water molecule in each state. The second set of simulations, referred to as MDQT-B, is carried out specifically to compare the time evolution of the vibrational state populations with previous theoretical results and to discuss the importance of the intermolecular VV energy transfer processes, and the quantum subsystem has in this case only two states, the ground vibrational state and the excited state corresponding to the water molecule whose bending mode is initially excited to the v = 1 level. For the sake of simplicity, we hereafter refer to the water molecule initially excited in the bend as the solute molecule and to the rest of the water molecules as the solvent. 2.2. Computational Details. We have used the flexible SPC/E model proposed by Ferguson47 to describe the H2O molecules. This model is extremely popular in molecular dynamics simulations because of its simplicity48 and has been shown to provide the best overall description of different vibrational spectroscopy experiments.49

The MDQT simulations were carried out by placing 256 molecules of water in the unit cell. The simulations therefore involved a number of 1536 classical translational and rotational coordinates and 256 and 1 quantum vibrational bending coordinates respectively for the MDQT-A and the MDQT-B sets of calculations. The length of the unit cell length was appropriately selected so as to reproduce the experimental water density at the desired temperatures (see Table 1) with periodic boundary conditions, and the Ewald sum method was used for the long-range interactions. Hamilton’s equations for the translational coordinates were integrated using the leapfrog algorithm,50 and the classical rotational motions were described using quaternions50 and propagated using the midstep implicit leapfrog rotational algorithm proposed by Svanberg.51 For every temperature, the system was first equilibrated during 300 ps, and then equilibration continued for a period of 1600 ps more, during which initial conditions were taken at time intervals of 2 ps. During equilibration, the temperature was maintained fixed at a given mean value by coupling to a thermal bath,52 and the quantum wave function was held constant at the vibrational ground state of the system. The MDQT simulations were carried out using a code developed by our research group. As shown below, different H-bonded structures can interconvert on the time scale of the vibrational relaxation. For this reason, we carried out some hybrid quantum/classical simulations keeping the time-dependent wave function fixed in the ground vibrational state in order to obtain information about the equilibrium structure of the water molecules and the H-bond network at different temperatures. These simulations were performed in the NVE ensemble in order to avoid any influence of the velocity scaling on the results. In the relaxation simulations, one random water molecule was instantaneously excited to its bend fundamental at t = 0, and the system was then propagated in the NVE ensemble. Since the energy transferred from the solute to the solvent (∼1400 cm−1) is much smaller than the total kinetic energy of the system, the temperature only increased in 2.6 K during the relaxation period. A classical integration step of 1.5 fs was used both in the whole equilibration process and in the relaxation simulations. The simulations were propagated only up to 1.2 ps during the relaxation process for computational demands. 9429

dx.doi.org/10.1021/jp5058447 | J. Phys. Chem. B 2014, 118, 9427−9437

The Journal of Physical Chemistry B

Article

3. NUMERICAL RESULTS 3.1. Equilibrium Dynamics. We have performed independent sets of simulations at the 277, 295, 323, and 348 K temperatures selected by Ashihara et al.11 in their experiments. In Table 1, we give the densities of the system at the different temperatures, which decrease as observed by 2.5%, as well as the values of the average number of H-bonds (nHB) and the number of molecules in the first (nfss) and in the second solvatation shells (nsss). Among the multiple ways of defining a H-bond,53 we use the simplest one which considers only the H···O intermolecular distance r. The first minimum of the radial distribution gOH(r) located at r = 2.41 Å is taken accordingly as the cutoff distance, so every intermolecular H···O pair is classified as H-bonded if r is less than the cutoff and as not H-bonded otherwise. We also assume that the radii of the first and second solvatation shells, hereafter referred to respectively as fss and sss, are those calculated by Voth et al.54 (rfss = 3.3 Å and rsss = 5.5 Å) using the SPC/E water model. At 277 K, the MDQT equilibrium simulations give a number of 3.83 H-bonds, which reduces to 3.63 at the highest temperature of 348 K, with a variation therefore of −5.2%, which is about twice that observed for the density. As for the number of molecules in the fss, there is a maximum at 295 K and a nonproportional variation, therefore, to the increasing temperature, which does not exceed 1.4%. In contrast, the number of molecules in the sss decreases progressively as the temperature rises (and the density gets lower), with an overall variation of −4.6%, similar to that observed for the nHB values. The number of H-bonds determined at 295 K, 3.78, is the same as that found for the HOD/D2O system32 and slightly higher than that obtained by Skinner et al.53 at 300 K, 3.7, who also used the SPC/E model. This small difference presumably arises from the use of a flexible model and from the quantum description made in this study of the vibrational bending motions of the water molecules. The numbers of molecules calculated in fss and sss at 295 K, 4.32 and 18.21, respectively, compare well, in turn, with those obtained by Voth et al.,54 4.34 and 18.25, at 298 K. The average number of H-bonds in the solute molecule provides no dynamical information about the structure of the liquid. To get that information, we use the instantaneous number of hydrogen bonds surrounding each water molecule, which we refer to as the H-bond environment. The relative importance of the H-bond environments and their lifetimes gives us, moreover, a deeper understanding of the molecular structure of the liquid. Table 1 additionally includes the percentage of water molecules with 1−2, 3, 4, and 5−6 H-bond environments. The 1 and 2, and the 5 and 6, H-bond environments are gathered together because they have small contributions which make their individual values difficult to converge. As observed, the H-bonds environment with the largest contributions corresponds to molecules with 3 and 4 Hbonds, versus the hypocoordinate 1−2 and hypercoordinate 5− 6 H-bond environments, which have smaller, although not negligible, contributions. The coexistence of different H-bond environments is supported by the experimental infrared absorption spectrum of H2O, which presents a bending absorption band with a pronounced shape anisotropy,14 attributed to two frequency contributions. From the variation of the bending band with the temperature, Walrafen et al.55 concluded that this anisotropy arises from the contribution of

two unresolved bending components: a high-frequency contribution due to the water molecules being fully H-bonded and another low-frequency broad component coming from the partially H-bonded molecules. We also see in Table 1 that as the temperature rises, the percentage of water molecules with 4 H-bonds decreases and the contributions of the remaining H-bond environments simultaneously increase. This is a consequence of the increase of the molecular kinetic energy, which allows for larger molecular displacements, favoring the appearance of the environments with 1−2, 3, or 5−6 H-bonds which are more energetic than the 4 H-bond environment. We also give in Table 1 the lifetimes of the different H-bond environments. The general trend is that the lifetimes of H-bond environments become shorter as temperature rises, but not at the same rate. Within the range of temperatures considered, we observe relative lifetime decreases of −42.7%, −32.5%, −24.2%, and −14.0% for the 4, 5−6, 3, and 1−2 H-bond environments, respectively. We have also estimated the average lifetime of the H-bonds environments using the percentages of H-bond environments and their individual lifetimes. These lifetimes decrease (see Table 1) from 78 to 47 fs as the temperature of the system rises from 277 to 348 K. We should emphasize that this decrease of the average lifetimes of the H-bond environments results from the combination of two effects: the lifetime decrease of all of the H-bond environments and the contribution increase of the 1−2 and 3 H-bond environments with respect to the 4 H-bond environment, which have the largest lifetimes at the temperatures considered. The calculated average H-bond environment lifetimes at 277 and 295 K, 78 and 67 fs, are in good agreement with the experimental value of 77 fs at 298 K given by Elsaesser et al.14 and also with the theoretical value of 80−90 fs at 300 K provided by Skinner et al.53 Table 1 finally includes the H-bond lifetimes, which are approximately 7 times higher than the lifetimes of the H-bond environments in liquid water at all the temperatures studied. We notice here that the lifetime of a H-bond environment is conditioned by the formation and breaking of any of the Hbonds around the water molecule and that there are however H-bonds that may endure for long times, so the average lifetime of an individual H-bond is expected to be much longer than the lifetime of the whole environment. We have also analyzed the frequency distribution of the vibrational bending excitations extracted from the MDQT simulations, which are calculated as the differences between the diabatic ground and excited bending energies. In Table 2 we show the mean bending excitation frequencies at different temperatures. At 295 K we see that the frequency distribution is centered at 1608.1 cm−1, a value which is ∼50 cm−1 lower than the experimental values measured at 298 K by different groups of 1644.5,56 1650,14 and 1645 cm−1.11 As for the width of the Table 2. Average Bending Excitation Frequencies (in cm−1) of the H-Bond Environments at Different Temperatures

9430

H-bond environment

277 K

295 K

323 K

348 K

2 3 4 5 mean value

1598.5 1606.9 1610.8 1611.6 1609.6

1598.1 1605.6 1609.4 1610.2 1608.1

1596.8 1604.8 1608.3 1609.2 1606.0

1596.9 1604.1 1607.9 1608.1 1605.1

dx.doi.org/10.1021/jp5058447 | J. Phys. Chem. B 2014, 118, 9427−9437

The Journal of Physical Chemistry B

Article

with the quantum population during the first 0.4 ps and decay slightly faster at longer times. The sum of the classical populations of the excited bending states of the solvent reaches a maximum of 40% at 300 fs and then decays approaching zero at longer times. A similar result is obtained at the other temperatures studied, which is indicative of the importance that the VV channel has in the MDQT-A simulations and of the independence of this pathway with the temperature. The classical population of the ground state grows continuously toward the unity value. The disagreement between the quantum and classical populations is due to the trajectories that become trapped after hopping to the ground state, since most of the transitions from the ground to the excited states are energetically forbidden.26,57 Consequently, only the classical population reproduces the Boltzmann thermodynamic behavior at the long time limit. Following previous studies by Thrular et al.58 on nonadiabatic electronic transitions, we have ignored the frustrated hops and left the classical momenta unchanged. The results extracted from the MDQT-B simulations show first that the population of the initial quantum state decays in this case much more slowly than in the MDQT-A simulations, which highlights the importance of the intermolecular couplings between vibrationally excited states. It is also observed that the classical population of the initial state diverges from the corresponding quantum population from the beginning. This is due to the fluctuations of the quantum population observed in the individual trajectories.31 These oscillations cause transitions to the ground state that are mostly irreversible, as already discussed above. The relaxation lifetimes, defined as the time required for the population of the initially bending state to decay by a factor of 1/e, are presented and compared with experimental values in Table 3. There are two fundamental differences between the

frequency distribution, the MDQT simulations provide values of 35.5 cm−1 at 277 K and 40 cm−1 at 295 K, narrower than the experimental bandwidths56 which vary from 88 cm−1 at 298 K to 83 cm−1 at 343 K due to the contribution of the homogeneous broadening. The rise of the temperature in the simulations from 277 to 348 K results in a mean frequency decrease of 4.5 cm−1, that is, in a band red-shift. We notice in this respect that Falk et al.56 have measured frequencies of 1644.5 and 1642.3 cm−1, at 298 and 343 K, respectively, with a decreasing ratio of 0.049 cm−1 K, which agrees well with that of 0.057 cm−1 K extracted from the MDQT simulations. Table 2 also includes the average bending excitation frequencies for the different H-bonds environments. As seen, the frequencies at every temperature increase with the number of H-bonds of the molecule, in good agreement with the blueshift of the experimental bending band observed from gas phase to liquid water.56 In addition, the frequencies of a given H-bond environment decrease as the temperature rises. We note, however, that the red-shift of the mean value of the band (−4.5 cm−1) is almost twice the red-shift observed for the frequency of every H-bond environment. This is due to the fact that the relative importance of the environments with fewer H-bonds, which have lower frequencies, increases as the temperature rises as previously discussed (see Table 1). The change of the distribution of the different H-bond environments accounts, therefore, for approximately half of the calculated red-shift of the band. 3.2. Vibrational Populations. In Figure 1 we show the time evolution of the populations of the vibrational states

Table 3. Lifetimes (in fs) of the Vibrational Relaxation of the Water Bend Fundamental temperature (K) MDQTa exptl

Figure 1. Time evolution of quantum (black lines) and classical (red lines) populations of the initial excited state (solid lines) and the ground state (dotted lines) obtained in (a) MDQT-A and (b) MDQTB simulations at 295 K. The dashed lines in (a) correspond to the sum of the populations of all excited vibrational states populated through VV intermolecular transfers.

A B Ashihara et al.b Elsaesser et al.c Lindner et al.d

277

295

323

348

214 574 176

220 552 174 170 220

224 498 200

231 507 250

Results with ±4 fs standard deviations. bReference 11: results with ±15 fs standard deviations. cReference 14: results with ±30 fs standard deviations, obtained at room temperature. dReference 17: results obtained at 25 °C.

a

involved in the relaxation of the H2O bending fundamental extracted from the MDQT-A and MDQT-B simulations at 295 K. The results from the simulations performed at the other temperatures probed, 277, 323, and 348 K, exhibit similar trends. Since the vibrational frequencies of the H2O bends are much higher than the thermal energy at room temperature, it is expected according to the Boltzmann distribution that only the ground vibrational state will be populated at equilibrium. As observed, in the MQDT-A simulations, the quantum population of the initial state decays steadily during relaxation, approaching zero at long times. Simultaneously, the sum of the quantum populations of the excited bending states of the solvent grow while the quantum population of the ground state always remains small with values around 0.2%. We see also that the classical populations of the initial excited state agree well

MDQT-A and MDQT-B lifetime values. The first is that the MDQT-B lifetimes are substantially longer, by a factor of about 21/2, than the MDQT-A lifetimes. The second difference is that the MDQT-A lifetimes increase with temperature, in contrast with MDQT-B lifetimes which decrease toward a constant value at the higher temperatures. The lifetime obtained in the MDQT-A simulations at 295 K of 220 fs is in reasonable agreement with the experimental values at similar temperatures, which vary between 170 and 220 fs,11,14,17 whereas the lifetimes form the MDQT-B simulations, in which the VV channel is closed, double these values. It is therefore clearly necessary to include the intermolecular VV channels in the simulations in order to obtain lifetimes consistent with the experimental data. The significant decrease of the relaxation lifetime when the VV 9431

dx.doi.org/10.1021/jp5058447 | J. Phys. Chem. B 2014, 118, 9427−9437

The Journal of Physical Chemistry B

Article

channel is open was previously stated by Hynes et al.18 in their classical simulations reporting a relaxation time of 150 fs at 300 K when including the VV channel versus the value of 270 fs obtained by freezing vibrational motion of the solvent. However, we should note that the contribution of the VV energy transfers in these classical simulations was not evaluated, so the modification of the relaxation time could be also related to the different couplings emerging from the rigid or flexible models of the solvent molecules. We also note that the increase of the MDQT-A lifetimes with temperature is similar to the increase observed in the classical simulations,25 which vary from 270 to 285 fs when the temperature increases from 300 to 360 K, while the experimental results by Ashihara et al.11 provide a faster increase. In our previous work on the relaxation of the bending water mode conducted using the Ehrenfest approximation,32 we obtained a relaxation time of 120 fs, about half that of the MDQT-A time obtained here. This is probably due to the introduction of quantum correction factors in the propagation of the time-dependent wave function32 needed to force the right thermodynamic behavior of the quantum populations at long times. Figure 1 shows also that the decay of the population of the initially excited states is not monoexponential, in agreement with our earlier results using the Ehrenfest method.20 In fact, the populations cannot be well reproduced using simple kinetic schemes. The reason for this seems to be the important role that the H-bond network surrounding the solute molecule plays during the relaxation process. We notice in this respect that time-resolved infrared pump−probe spectroscopy measurements of the OH-stretch vibrational relaxation of HOD in liquid D2O at temperatures up to 700 K carried out by Schwarzer et al.5 concluded that the time scales of spectral diffusion and vibrational relaxation in liquid water at room temperature are apparently not well separated, with both processes occurring in parallel. If the decay of the initial bending state population depends on the number of H-bonds surrounding the excited molecule due to an alteration in the potential couplings that induce the relaxation, different lifetimes should be obtained in the MDQT simulation depending on the initial number of H-bonds. In Table 4, we give the relaxation

relaxation process, as previously proposed by Hynes et al.18 and supported by the experimental data.59 We also observe in Table 4 that the relaxation lifetimes for a given initial H-bond environment tend in general to increase with the temperature, although there are some fluctuations. Interestingly, the increases observed for every H-bond environment are smoother than those found for the mean values of the relaxation lifetimes. The reason for this is the second contribution that arises from the change in the relative importance of the different H-bond environments with the temperature (see Table 1). We have also found that excitation of molecules with fundamental bending frequency below the mean values (see Table 2) gives slower relaxation lifetimes than those of the molecules whose excitation frequencies lie in the blue-shift of the distribution. This is due to the relationship between the number of H-bonded molecules and the transition frequency shown in Table 2. Molecules with lower transition frequencies are to be associated with hypocoordinated environments and therefore exhibit slower relaxation lifetimes. We should emphasize that the lifetimes included in Table 4 cannot be considered, in general, as those corresponding to a single particular H-bond environment. As observed in Tables 1 and 4, the lifetimes of the H-bond environments are significantly shorter than the lifetimes of the vibrational relaxation. The molecule initially excited therefore modifies its H-bond environment several times during the vibrational relaxation, thus smoothing the differences of the lifetimes shown in Table 4. 3.3. Quantum Hops and Vibrational Energy Redistribution. The capacity of the MDQT method to give detailed information on the vibrational hops during the relaxation dynamics allows us to quantify the relative contributions of the H2O inter- and intramolecular channels to the bend relaxation. We consider the seven types of hops occurring during the relaxation process previously proposed for the HOD/D2O system.32 They are the types 1, 3, and 5 of hops, which account respectively for the relaxation of the initially excited H2O molecule to the ground state, the intermolecular VV energy transfer from the initially excited H2O molecule to another H2O solvent molecule, and the relaxation of the solvent H2O molecules to their ground vibrational states. Next we have hops 2, 4, and 6 corresponding to the inverse processes and finally the hops of type 7 accounting for the VV energy transfer between two solvent H2O molecules (see eqs 14−17 in ref 32). From the MDQT simulations we obtain then that hops of type 2 have, as expected, a very small contribution (