Molecular dynamics calculation of the infrared spectra in liquid H2O-D2O mixtures

Journal of MolecularLiquids,62(1994) 17-31

MOLECULAR DYNAMICS CALCULATION OF THE INFRARED SPECTRA IN LIQUID H20-D20 MIXTURES J.Martl(a'b), J.A.Padr6(a) and E.Guhrdia(b) (a)

(b)

Departament de Flsica Fonamental. Universitat de Barcelona. Diagonal, 647. 08028 Barcelona. Spain. Departament de FIsica i Enginyeria Nuclear. Universitat Politecnica de Catalunya. Sor Eulalia d#Anzizu, 08034 Barcelona. Spain. (Received 20 July 1993; accepted 28 February

B4-B5.

1994)

Abstract the distributions corresponding to The frequency vibrational and librational molecular motions in liquid water isotopically substituted mixtures are investigated. The dynamics spectra is calculated from molecular infrared simulations assuming a flexible single point charge (SPC) potential. A new set of parameters for the intramolecular

and

potential which allows a reliable reproduction of the peaks of the experimental infrared spectra in the gas and liquid phases is proposed. Pure H20, D20 and HDO liquids as well as dilute HDO solutions are analyzed. It has been verified that the OH and OD stretching modes in HDO are very close to those in H20 and D20 and that the HDO spectra is little influenced by different environments (HDO, H20 and D20).

1. INTRODUCTION Molecular dynamics (MD) simulation provides detailed pictures of the microscopic motions in liquids. This kind of information can be used for the study of the spectroscopic

‘0167-7322/94/$07.00 SSDI 0167-7322

0 1994 - Elsevier Science B.V. All rights reserved. (94) 00769-l

18

properties of molecular liquids [l-4]. Though it is a classical approximation, MD is a helpful tool for the interpretation of the the results of spectroscopic experiments. Moreover, comparison of MD results with experimental data can be used for testing the potential models employed in the simulations. The experimental determination of the infrared (IR) spectra of liquid H20 and D20 is a difficult task due to the large absorption coefficients of these substances. Moreover, and the these spectra are broad, overlap considerably contributions of the different vibrational modes cannot be easily determined because of the coupling effects (Fermi resonance) between the fundamental stretching modes and the first overtone of bending [5,6]. A useful alternative is the study of the IR spectra of HDO in dilute mixtures of H20 in D20 or

vice

versa.

The main advantage of considering the HDO molecules is that the three fundamental modes and many of the overtones are widely separated [6,7]. Several experimental IR studies of the vibrational spectra of dilute solutions of HDO in H20 and D20 were carried out [5-101. However, little attention has been paid to the MD study of the HDO spectra. It should be noted that, unlike in real experiments, we can perform MD simulations of pure HDO which is equivalent to a dilute solution of HDO in HDO. This allows us a substantial improvement in the statistic of the simulation results. MD calculations of the IR spectra of H20 and D20 using different potentials have been recently reported [3,11,12]. In this paper, we apply the MD simulation technique to the analysis of the translational and librational frequency bands in H20, D20 and HDO. Previously we have refined the parameters of the intramolecular potential Toukan and Rahman [4] so that the infrared from MD become closer to the experimental differences between the vibrational bands

model proposed by spectra resulting observations. The obtained for the

three systems are analysed and the influence of the environment on the IR spectra of HDO is investigated.

19

2.- COMPUTATIONAL

DETAILS

2.1 - Interaction

Potentials

The potentials used in the MD simulations were the SPC for the intermolecular model of Berendsen et al. 1131 interactions and the potential by Toukan and Rahman [43 with the form detailed in for the intramolecular ref.14 interactions. Previous MD simulations [3,4,15] of H20 and D20 using the Toukan and Rahman intramolecular potential showed that the positions of the bending peaks differ from the experimental values about 150-200 cm-' (in the case of the stretching the accordance was quite good). In order to achieve a better reproduction of the experimental data we have used a trial and error procedure to find a new set of parameters for the

intramolecular potential. The resulting parameters are similar to those in ref.[4] (see Table l), but, as will be shown in Section 4, the frequency peaks resulting with these new parameters are notoriously closer to the experimental ones for both bending and stretching in liquid and gas phases.

Table Parameters

1

of the intramolecular

Ref. [4] ke (mdyn A-')

I

’

potential

*

This work

2.283

1.697

kr,9 (mdyn A-')

- 1.469

- 1.280

krrr(mdyn A-')

0.776

0.880

D (mdyn A)

0.708

0.667

a (K-l)

2.567

2.567

The potential form is given in the Appendix of ref.[l4]

20

2.2.- MD Simulations MD simulations of pure H20, D20 and HDO as Well as of HDO-H-0 and HDO-D-0 mixtures were carried out. We considered 216 mtlecules (in the case of mixtures 16 HDO + 200 H20/D20) in a cubic box with the usual periodic boundary conditions. The size of the box was chosen so that the density was p = 1 g/cm3 and the temperature was kept at T = 298 K. A leapfrog Verlet algorithm with coupling to a thermal bath [16] with an integration time-step of 0.5 fs was employed. The equilibration of the translational and intramolecular motions were performed separately [17]. The flexible SPC potential model of Toukan and Rahman [4] with the two set of parameters reported in Table 1 was assumed. The Ewald method was used for the long-range coulombian interactions. In order to calculate the IR spectra of pure water corresponding to the gas phase we carried out short MD runs of a single molecule with a kinetic energy corresponding

to

the

room

These fast MD temperature. simulations of water monomers were also very useful for fitting the intramolecular potential parameters. 2.3.- Calculation of the IR Spectra Two procedures were used for the calculation of the infrared spectra: i) The most common method is based on the determination of the velocity autocorrelation functions of the atomic velocities (CH(t)I CD(t) I Co(t)) during the MD simulations. Then, the IR spectra (G:(U), G:(U), S:(U)) are calculated as the Fourier transform of the corresponding velocity autocorrelation functions. ii) The second method [3,12] is based on the decomposition of the instantaneous velocities of the H or D atoms relative to the molecular center of mass (ui(t)) into components parallel and perpendicular to the molecular plane ui(t) = u!(t) + I_+)

;

i = 1,2

where i refers to the hydrogen or deuterium atom of a

(1) molecule

21

FIGURE l.- Unitary vectors and rotation axes used in this work.

and

u;(t)

=

u;(t)+)

(see

Fig.1).

Moreover,

u:(t)

is

decomposed into components parallel (Bi(t) = Bi(t)ei(t)) and perpendicular (Pi(t)= Pi(t)fi(t)) to the OH or OD directions . Then, three combinations of the resulting velocity components are defined. In the case of Hz0 or D20 Ql(t)=BI(t)+B2(t); QZ(t)=Pl(t)+P2(t); Q3(t)=Bl(t)-B2(t)

(2)

Ql(t) and Q,(t) are combinations of atomic motions along the OH directions which approximately describe the symmetric and asymmetric stretching modes, respectively. Q,(t) is a combination of velocities which are perpendicular to OH or OD and it corresponds to the bending mode. In the case of HDO Q,(t) = BI(t); Q,(t) = P,(t)+P,(t); Q,(t) = B2(t) where Q,(t) and Q,(t) describe the stretching modes, respectively, and

hydrogen Q2

and

corresponds

(3) deuterium to

the

bending.

Although rotations are strictly defined for rigid bodies, it is usual to describe the rotational motions of flexible molecules as a combination of rotations around characteristic molecular axes. The procedure proposed for molecular vibrations

22 may be extended to rotations [3,12] following time dependent quantities R,&t)=Pl(t)-P2(t);

Ry(t)=u;(t)+u;(t);

by

considering

RZ(t)=u;(t)-u;(t)

the

(4)

where the X-axis is perpendicular to the instantaneous plane of molecule and Ry(t) and RZ(t) are associated with the rotations around axes in the molecular plane (2 is along the dipole moment)

(see Fig.1).

The

IR

motions

spectra

are

corresponding

calculated

by

autocorrelation

functions

defined

(2), (3) and (4).

by Eqs.

3.- RESULTS

of

the

to

the

Fourier Q(t)'s

different

kind

transforming and

R(t)'s

of the

variables

AND DISCUSSION

3.1.- Gas Phase: H20 The vibrational intramolecular motions of water in the gas phase have been analysed from MD simulations of a single

molecule using both the original potential parameters of Toukan and Rahman and those proposed in this work (Table 1). We have determined the stretching and bending spectra according to the second method described in Sect.2.3. The results are summarized in Table

2 and they show that when the new parameters

are

employed the agreement among the frequencies of the vibrational peaks from MD and from experiments [18] is quite good. We have also calculated the maxima of the frequency spectra for the two sets of potential parameters using the harmonic analysis method [19]. The results are shown in Table 2. It should be noticed that the harmonic frequencies for the stretching modes are in accordance with the frequencies of free water monomers. However, marked discrepancies are found in the case of bending. This indicates that the anharmonicity effects associated to the potential model are noticeable even in the gas phase. This may be the reason for the discrepancies between

23

the bending frequencies of liquid H20 and D20

[3] resulting

from MD using the potential parameters of Toukan and Rahman [4] and the experimental data (it should be remembered that the parameters in [4] were obtained from the harmonic frequencies). Table 2 -1 Vibrational frequencies (cm ) (Gas Phase)

:P1::5E

-;::,-

rhis work r-R[4]* Bending

1594

1520

1725

1433

1652

Sym.Stretch.

3656

3685

3815

3683

3820

Asym.Stretch.

3755

3785

3945

3798

* Using the potential parameters of ref [4]

3.2.- Pure Liquids: H20, D20 and HDO The spectra of the liquid phases have been calculated using the two methods described in Sect. 2.3. The resulting S!(U) and it(~) are represented in Fig.2. The librational and stretching bands of H and D in HDO are very close to the corresponding bands in H20 and D20, respectively. This corroborates the validity of the assumption in experimental works that the IR stretching frequencies from HDO can be extended to the H20 and D20 molecules. As expected, in the case of bending, the results for HDO are markedly different from those for symmetric H20 or D20 molecules and it may be observed that the bending frequencies decrease when the molecular masses increase. The symmetric and asymmetric contributions to the stretching bands of H20 and D20 (Fig.3) have been determined

24

(a)

FIGURE

2.-

Spectra

of

o/10'

(cm-‘)

o/10'

(cm-‘)

the

hydrogen

and

deuterium

velocity

autocorrelation functions. a) hydrogen in H20 ( ) and HDO ) and HDO (------). (------) b) deuterium in D20 (

from the time autocorrelation functions Q,(t) and Q,(t) defined in (2). In the case of HDO, the spectra obtained from Q,(t) and Q,(t) correspond to the stretching bands resulting from CH(t) and CD(t), respectively. The i:(o) spectra calculated from the simulations of H20, D20 and HDO are represented in Fig.4. The three functions show a noticeable maximum around 50 cm-l and a shoulder around 200 cm-1 which have been associated to the motions between bending stretching intermolecular and hydrogen-bonded molecules, respectively [lo]. It should be noted that the shoulder becomes less pronounced when the molecular mass increases.

25

0.8

D20

HDO 0.8

0.0 2.0

2.5

3.0 o/10’

FIGURE

3 .-Spectra

3.5

4

(cm-')

of Q,(t) and Q,(t) as defined

in (2) and

(3).

a) Symmetric

(-

) and asymmetric

(------) stretching

in H20

b) Symmetric

(-

) and asymmetric

(------) stretching

in D20

c) Hydrogen

and deuterium

stretching

in HDO.

26

w/10'(cm-')

FIGURE

4.-

Spectra of functions in H20 (-),

the oxygen velocity autocorrelation D20 (------) and HDO (- - - -).

Table 3 -1 Vibrational frequencies (cm ) (Liquid Ii20and DZO) Heavy Water

Water

Experiment* MD Simulat. Experiment* MD Simulat, Bending

1645 (100)

1650 (280)

1210 (60)

Stretching

3390 (375)

3385 (455)

2500 (300) 2470 (335)

Sym.Stretch.

---

3370 (440)

_--

i\sym.Stretch.

---

3410 (355)

---

*

1200 (180)

2365 (285) 12505 (255)

From Ref. [20]. Values in parenthesis are the band widths at l/2 maximum.

27

Since the spectral bands are rather broad the positions of the peaks are not well defined. However, we have determined the frequency peaks for the different vibrational modes with an -1 and they are summarized in Tables 2 accuracy of about f10 cm and 3. The results for H20 and D20 [20] are in very good agreement with ,the experimental data. This confirms the suitability of the new parameters for the intramolecular potential. Moreover, the accordance between MD results in both gas and liquid phases shows that the shifts associated to the intermolecular interactions in liquid water and heavy water are realistically reproduced by the assumed SPC potential. In the case of HDO, the stretching bands are slightly closer to the experimental findings [9] when the original parameters are used but there is a clear improvement in the bending band when the new parameters are considered. It should also be noticed that as the librational bands are mainly influenced by the intermolecular interactions, they are not affected by the small changes in the intramolecular motions.

Table 4 Vibrational frequencies (cm-1 ) (Liquid IWO)

T Bending

*

MD Simulation Experiment*

This work

1450 (85)

1450 (195)

H-Stretching

3440 (250)

3355 (345)

D-Stretching

2530 (170)

2440 (250)

1620 (150)

From ref. [9]; ** Using the potential parameters of ref.[4]. Values in parenthesis are the band widths at l/2 maximum.

28 We have also compared the widths of the frequency bands at l/2 of the maxima of our spectra with those from the experiments (Tables 3 and 4). The agreement is only qualitative and the MD bands are wider in all cases. This should be mainly associated to the use of a classical description.

3.3.- Isotopic mixtures: HDO in Ii20and D20 In contrast to the MD findings presented in Section 3.2, the experimental IR spectra of HDO were obtained from dilute solutions of HDO in H20 and D20 [7-lo]. In order to analyze the influence of the environment on the HDO spectra we also carried out MD simulations of dilute solutions of HDO in H20 and D20. The results of these simulations corroborate that the HDO spectra do not show significant dependence on the considered environments. The bending and stretching frequency bands are practically identical. In the case of libration, which is an intermolecular band, we have observed small differences. In order to investigate these differences we have decomposed the librational band into components corresponding to the rotations around the three molecular axes defined in Section 2.3 (see Fig.1). As may be observed in Fig.5 the biggest differences among the rotational spectra for the three mixtures correspond to the Y axis. This shows that the influence of the environment is bigger for the rotation corresponding to a smaller moment of inertia of the HDO molecule.

4 . - CONCLUDING REMARKS The results of this work show that classical MD simulations using a suitable potential can acceptably reproduce the positions of the peaks in the IR spectra of water and isotopically substituted water mixtures. When the SPC potential [4] with the intramolecular parameters proposed in this paper is used, both the MD stretching and bending frequencies for liquid H20 and D20 are in good agreement with the experimental

29

X-axis

0.8

0.0

0.3

0.8

0.5

1

w/10' (cm-')

Y-axis '\ 0.8

0.0

0.3

0.5

0.8

1

o/10' (cm-')

Z-axis

0.0

0.3

0.8

0.5

o/ 10’

(cm-‘)

FIGURE 5.- Spectra of RX(t), Ry(t) and Rz(t), as defined in HDO in Hz0 (------) and HDO in D20 Eq.(4). Pure HDO (-), (- - - -).

30

data. In the case of the asymmetric HDO molecules the agreement is also quite good. Moreover, the frequency differences between the gas and liquid vibrational spectra of water are also well reproduced. It has also been shown that the stretching frequencies at a higher temperature (523K) are in accordance with experiments [21]. However, marked disagreements have been found in the width of the frequency bands. It should be noted and comparison of the MD simulation that a complete experimental findings cannot be based on the spectral densities obtained in this work but the imaginary part of the dielectric constant should be calculated [22,23]. Nevertheless, a detailed reproduction of the shape and intensity of the experimental infrared spectra would require to perform quantum mechanical simulations. We want to emphasize the usefulness of the MD simulation technique for testing different approximations. It has been corroborated that the OH and OD stretching modes in liquid HDO do not show significant differences with those in liquid H20 and D20. Moreover, it has been verified the HDO spectra is little influenced by different environments (HDO, H20 or D20). Finally, it should be pointed out, that MD simulations provide reliable decompositions of the libration and stretching bands which cannot be rigorously deduced from the experimental data.

ACKNOWLEDGMENTS work has been supported by "Centre de This the Supercomputaci6 de Catalunya" under a CESCA grant. Financial support of DGICYT, project PB90-0613-CO3 is also acknowledged.

31

REFERENCES [l] Madden,P.A. in Molecular Dynamics Simulation of Statistical Mechanical Systems, 1986, Ed. by G.Ciccotti and W.G.Hoover (North- Holland). [2] Evans,M.W., 1992, Adv.Chem.Phys. al, 361. [3] J.MartfL, J.A.Padr6 and E.Gu&rdia, Mol.Simul. (in press). [4] Toukan,K. and Rahman,A. 1985, Phys.Rev.B, 31, 2643. Conway,B.E., 1981, Ionic Hydration in Chemistry and r51 Biophysics, (Elsevier), Ch.7. and 161 Eisenberg,D. and Kauzmann, 1969, The Structure Properties of Water, (Clarendon), Ch.4. [7] Falk,M. and Ford,T.A., 1966, Can.J.Chem. 44, 1699. [8] Waldron,R.D., 1957, J.Chem.Phys., 26, 809. [9] Wyss,H.R. and Falk,M., 1970, Can.J.Chem., 48, 607. [lo] Walrafen,G.E., 1971, In Water, A Comprehensive Treatise, Franks,F. Ed., Vol. 1 (Plenum), Ch.5. [ll] Bansil,R., Berger,T., Toukan,K., Ricci,M.A. and Chen,S.H., 1986, Chem.Phys.Lett., 132, 165. [12] PZilinkas,G., Bak6,1., Heinzinger,K. and Bopp,P., 1991, Molec. Phys. 73, 897. [13] Berendsen,H.J.C., Postma,J.P.M.,van Gunsteren,W.F.,Hermans, J. in Intermolecular Forces, 1981, Ed. by B.Pullman (Reidel). [14] Gu&rdia,E. and Pad&,J.A., 1990, Chem.Phys. 144, 353. [15] Dang,L.X. and Pettitt,B.M., 1987, J.Phys.Chem. 91, 3349. [16] Berendsen,H.J.C.,

Postma,J.P.M., van Gunsteren,W.F.,

di

Nola,A. and Haak,J.R., 1984, J.Chem.Phys. al, 3684. [17] Wallquist,A. and Teleman,O., 1991, Molec.Phys. 74, 515. [18] Chen,S.H., Toukan,K., Chung,K.L., Price,D.L. and Teixeira, J. 1984, Phys.Rev. Lett. 53, 1360. [19] Wilson,E.B., Decius,J.C. and Cross,P.C., 1955, Molecular Vibrations (McGraw-Hill, New York). [20] Bertie,J.E., Ahmed,M.K. and Eysel,H., 1989, J.Phys.Chem. 93, 2210. [21] J.A.Padre, J.Marti and E.Guardia, 1994, J.Phys.Conden. Matter 6 (in press). [22] D.A.McQuarrie,Statistical Mechanics,l976,Haper & Row,Ch.21 [23] I.Ruff and D.J.Diestler, 1990, J.Chem.Phys. 93, 2032.