A molecular dynamics study of the OH stretching vibrational spectrum of liquid water

CHEMICAL PHYSICS LETTERS

Volume 132, number 2

A MOLECULAR DYNAMICS STUDY OF THE OH STRETCHING VIBRATIONAL R. BANSIL,

OF LIQUID WATER

T. BERGER

Center for Polymer Studies and Department

K. TOUKAN

SPECTRUM

5 December 1986

‘, M.A. RICCI

Nuclear Engineering

Department,

MIT,

of Physics, Boston University, Boston, MA 02215, USA

’ and S.H. CHEN Cambridge,

MA 02139, USA

Received 27 December 1985; in final form 16 September 1986 A classical molecular dynamics (CMD) study has been carried out to investigate the OH stretching spectrum of liquid water. The potential used is the simple point charge model modified to include anharmonic vibrational potentials. A series of CMD runs were made for a system of 100 water molecules. Spectral functions for infrared, Raman and inelastic neutron scattering were directly computed from the trajectories of the CMD. The computed spectral band shapes are in good qualitative agreement with experimentally measured spectra.

1. Introduction The OH stretching region of the vibrational spectrum of liquid water has been extensively investigated by both Raman and infrared [ 1,2] spectroscopy. Recently the same spectral region has been studied using incoherent inelastic neutron scattering spectroscopy [3]. The observed Raman and infrared band shape has traditionally been interpreted ln terms of some phenomenological models [ 1,2]. In recent years, a number of approaches have been developed to calculate the Raman and infrared stretching spectrum of liquid water based on rigid molecule classical molecular dynamics (CMD). By combining the CMD trajectories of rigid water molecules interacting via the ST-2 potential and normal-mode analysis Rice and co-workers [4,5] have calculated the infrared and Raman OH stretching spectra. They use the ST-2 potential to generate the configurations of 216 water molecules and then use either the instantaneous or time-averaged configurations as the input for a normal-mode analysis using an empirical vibrational potential which includes anharrnonicity and various coupling terms. An alternative approach has been developed by Watts and co-workers [6], where they used Monte Carlo calculations to generate the configurations of rigid water molecules ln the liquid state using the RWK2 potential and then performed a quantum local-mode analysis for a single water molecule using an intramolecular potential consisting of a sum of three Morse oscillators to obtain the infrared spectrum. A third approach is that of Postma et al. [7], who also used a rigid molecule CMD to first generate the configurations from which they calculated the instantaneous forces on a particular molecule due to its neighbors. This gave them a measure of the degree of distortion of the internal coordinates that the water molecule would have undergone due to the momentary environment of the molecule. By averaging this process over a sufficiently large number of configurations, they obtained the vibrational band shapes. An alternative approach to using rigid molecule CMD and some form of vibrational-mode analysis is to generate the trajectory directly with a flexible molecular potential [8] * which includes both inter- and intra-molecular ’ Permanent address: Research Institute, University of Petroleum and Minerals, Dhahran, Saudi Arabia. 2 Permanent address: Dipartimento di Fisica, Universitg di Roma “La Sapien&, Piazzale Aldo Moro, 00185 Rome, Italy. * For a modification of this approach see ref. [9].

0 009-26141861% 03.50 0 Elsevier Science Publishers (North-Holland Physics Publishing Division)

B.V.

165

Volume 132, number

2

CHEMICAL

PHYSICS

5 December

LETTERS

1986

terms. In this type of calculation the various vibrational spectra are obtained by Fourier transforming the appropriate time correlation functions. An obvious advantage of a flexible molecule CMD is that the vibrational motion occurs simultaneously with the translational and rotational motions and thus the coupling of these different types of motion is correctly accounted for. Lemberg and Stillinger [8] developed a central-force model that enabled both the intra- and inter-molecular potential to be described in a simple set of atom-atom interactions. They calculated the density of vibrational states for the proton motions and the infrared spectrum from the dipole-dipole correlation function. However, this approach was not pursued much further, because the central-force model gave a rather poor agreement with gas-phase properties and predicted vibrational frequencies were way off. Very recently Toukan and Rahman [lo] have rejuvenated this approach by using anharmonic Morse oscillators for the intramolecular part (as in ref. [6] ) and the SPC model [ 1 I] for the intermolecular part of the potential. They calculated the proton velocity autocorrelation function and its power spectrum and obtained much better agreement with experiment than in the previous calculations based on central-force potentials. In this Letter we discuss the implications of using a flexible anharmonic potential to see if the inclusion of anharmonicity produces better agreement for infrared and Raman spectra with experiment as well than calculations based on central-force potentials. In order to obtain a more comprehensive test of this model we calculate all three vibrational spectra: the Qdependent inelastic neutron spectrum (INS), the infrared (IR) spectrum and the Raman spectrum. These are obtained by Fourier transforming the proton current correlation function, the OH-bond dipole-dipole correlation function and the polarizability derivative correlation function. Strictly speaking, the calculation of dipole moments and polarizabilities involves a quantum calculation of the electronic states and the response of these states to an external electric field. Work along these ab initio lines has recently been possible for water monomers [12]. However, such a calculation for liquid water is at present not feasible. To get by this difficulty we have taken the same route as Belch and Rice [5] in using an empirical bond polarizability and bond dipole model to write the polarizability derivative tensor and the dipole moment in terms of atomic coordinates calculated in the CMD trajectory at each time step. Thus our approach combines certain aspects of the ReimersWatts infrared calculations [6] (the use of anharmonic Morse potentials) with those from the Belch-Rice cahlations [5] (the use of empirical bond dipole and polarizability models) with a flexible molecule CMD. Finally it should be noted that in the past all three vibrational spectra have never been calculated from the same model. However, by calculating all three vibrational spectra from the same model one can obtain a much better test of the validity of the model. We obtain very good qualitative agreement with experiment suggesting that it is crucial to include anharmonicity in the intramolecular potential to describe the vibrational spectra of liquid water. The calculations also indicate that in order to obtain better quantitative agreement one would have to include higherorder terms in the potential such as dipole-dipole coupling and Fermi resonance.

2. CMD method The intermolecular potential model employed in this CMD is the simple point charge model originally proposed by Berendsen et al. [ 111. The potential consists of a pairwise Lennard-Jones interaction between the oxygen sites, voO = -(A/@ t @/r)12 with.4 = 0.37122 nm (kJ mol-1)u6 and B = 0.3428 nm (kJ mo1-1)1112, together with Coulomb interaction between suitable point charges situated at the oxygen and hydrogen sites. We assign a charge +0.4 1 Ie I to each hydrogen site and -0.82 Ie I to the oxygen site. The intramolecular potential is an anharmonic Morse function for the O-H bond and a harmonic function for the H-H interaction. Internal coupling between the two O-H bonds, as well as coupling between the stretching and bending modes, has also been inserted in the potential: V mol ‘&H Here b1, 166

{[I

-

eXp(--p&l)]2

+ [1 -

eXP(-@2)]

2, + 4 b.@

+ c(LL’1 + b,)

b3

+ db,

b2.

(1)

Ar2 are the stretch in the OH bond lengths, and Ar3 is the stretch in the H-H distance. The parameters

Volume 132, number 2

CHEMICAL PHYSICS LETTERS

5 December 1986

DoH, p, b, c and d are taken from Toukan and Rahman ** WithDOH = 0.708 mdyn/A, p = 2.37 A-‘, b = 2.283 mdyn/A,c = -1.469 mdyn/A and d ~0.776 mdyn/A. These parameters reproduce the experimentally observed frequencies for a free water molecule and the dissociation energy for an OH bond. The model has recently been used to simulate water by Toukan and Rahman [lo] at three temperatures. They calculated the structural properties of water, which were similar to those previously generated by the “ST-Z” model [13]. In fact the partial pair correlation function,gHH(r), generated by the “SPC-flexible water” model was shown to be even closer to the neutron experiment [14] than the corresponding quantity from the “ST-2” model. They also calculated the velocity autocorrelation function (VAF) of protons and oxygens. It is well known that the Fourier transform of the VAF of protons is the proton density of states,fp(w), which is probably the most fundamental quantity describing the dynamics of hydrogen-bonded systems, such as water. Thefp(w) they generated indeed showed features similar to those of the experimentally measured INS spectra. The system we studied consists of 100 molecules in a constant-volume simulation cell with periodic boundary conditions. CMD runs at a density of 1 g/cm3 and temperature T = 250 K have been done with a time increment At = 6.5 X lo-l6 s. The system was aged for several thousand time steps to ensure reaching thermal equilibrium and the pair correlation functions were checked to see that the correct static structure was generated in the simulation [lo]. After equilibration the simulation was run for 9000 time steps, i.e. 5.8 ps. 2.1. Simuhtion of inekzsric neutron spectra Due to the large incoherent scattering cross section of hydrogen atoms as compared to oxygen, the INS primarily measures the proton dynamics. We shall focus our attention on the following single particle current density autocorrelation function, C&Q, t), defined as

, 7

where Uj and rj refer to the velocity and position of thejth proton in the system. The Q vector in eq. (2) is sampled from the relation: Q = (2n/L)@, v, 6), with (cl, v, 6 = 0, +l, *2 , ...). where L is the length of the cubic simulation cell. The choice of the Q values as given above is a consequence of the periodic boundary conditions imposed on the simulation cell. The scalar quantity C&Q, t) in eq. (2) is obtained by averaging over all possible orientations of the momentum transfer, Q, for a given Q value. The Fourier transform of the current autocorrelation function, G,(Q, w), can be shown to be related to the proton self-dynamic structure factor, S&Q, w), by [ 151

G&Q, a)= j- C&Q, T) cos(oOdt = (02/Q2)S,(Q,

0).

From a practical point of view G&Q, o) is the most natural quantity to calculate for the high-fre uency 1 modes of proton motions, because the low-frequency part of S,(Q, o) is attenuated by the factor o /Q2. This leads to’a rapid decay of the time correlation function C&Q, t), which allows the Fourier transform of this function to be computed from relatively short time histories. Furthermore, since lim,,,C,(Q, t) = (u&O) * uH(t)) and lime+OCs~Q, w) = fJw), the quantities C,@, t) and G&Q, w ) are natural finite-Q extensions of the proton velocity autocorrelation and the frequency distribution functions. Several preliminary test runs of CMD showed

** All the parameters except p are the same as In ref. [lo]. P was changed from 2.566 to 2.37 A +, because this accounts better for the damping in the OH Morse oscillator.

167

CHEMICAL PHYSICS LETTERS

Volume 132, number 2

5 December 1986

25

0.000 0

200

400

EC0

t/At

0

ZOO f/D

4ca

6Oa

0

2000 f&T-:,

4000

6000

(d) SCAllERlNGDEW-I? SPECTRUM Q=8.72 A-’

0 I!!L 0

*lo-'

2000

4000

Em0

WC-‘:,

Fig. 1. (a) Proton velocity autocorrelation function (VAF) at 250 K. Note that the time step interval & = 6.25 X lo-r6 s. (b) Proton current density autocorrelation function Cs(Q, r) calculated at 250 K for Q = 8.72 A-‘. (c) Proton density of states calculated by Fourier transforming the VAF. (d) Proton density scattering spectrum G,(Q, w) at Q = 8.72 A-‘. Spectral intensity attenuation and band broadening at the finite Q value is clearly seen from a comparison of (c) and (d).

that the current autocorrelation function C,(Q, t) decays to 1% of its initial value in less than 600 at. An 8000 time step trajectory was divided into four sets of 2000 time steps each to obtain a good average for the time correlation function at Q = 0 and Q = 8.72 A-‘. Time correlations were calculated over 1200 time step periods. The time origin was shifted by 20 steps and this shifting of the origin was done 40 times within each set. The time correlations from each of the four sets were then averaged to produce the final autocorrelation functions shown in figs. la and lb. Comparison between the two functions shows the quicker damping of the correlation function at the large Q value as mentioned previously. The frequency spectra for each of the two Q points were also calculated as shown in figs. Ic and Id. Comparing the two spectra, one notes the Q-dependent broadening of each of the three bands for large Q values as compared with the well-defined bands at the Q = 0 limit. In particular, this broadening is most pronounced in the librational band. Another aspect is the lower intensity of the entire spectrum at the higher Q value due to the damping at large momentum transfer. Furthermore in the finite-e spectrum, two features appear on the Stokes side of the bending and stretching bands at 2440 and 3760 cm-‘. These features look similar to the “satellite peaks” observed in the experimental neutron scattering data [3,16]. The same features have been reported recently in an atomic simulation by Lie and Clementi [ 171. A more detailed analysis of these weak bands will be discussed in a forthcoming work. 2.2. Simulation of IR spectrum The IR spectrum can be calculated from the dipole-dipole time correlation function, CIR(t) [18]. To evaluate this quantity we use the same model for the bond dipole moment as was used by Belch and Rice [5]. Since we are interested only in the time variation of the correlation function, we retain only the first-order time-varying term in the bond dipole model and thereby obtain,

(4) where N denotes the number of OH bonds, Gr$t) is a vector directed along the jth OH bond and of magnitude equal to the change in the bond length at time t from the equilibrium bond length. The dipole derivative, (a~/ari)o, is assigned the experimental value of 3.75 D/A [19]. To calculate the IR time correlation function we divided the 90OO:step CMD trajectory into three sets of 3000 steps each. In each set the IR time correlation function was evaluated over 1800 time steps, and was averag168

Volume 132, number 2

5 December 1986

CHEMICAL PHYSICS LETTERS

3. Discussion In fig. 2 we plot the infrared and the isotropic and anisotropic Raman time correlation functions as given by eqs. (4), (6) and (7) respectively. It is to be noted that the INS correlation function decays much faster than the Raman and IR correlation functions. Another obvious difference is the damped-oscillator-like behaviour of the optical cor-

(c) lSO-~~-uAMAN

CORRELATION (a)

(b)

10

,

‘“1*10”,

1

FUNCTIONS

INFRARED

ISOTROPIC

2800

3200

2m

3200

3600

4ow

3600

4000

RAMAN

I

03 2800

3200

3600 cm-’

4000

cm-’

,o e

(c)

ANISOTROPIC

RAMAN

Fig. 2. Correlation functions for (a) infrared and (b) isotropic and (c) anisotropic Raman calculated at 250 K using eqs. (4), (6) and (7), respectively. The correlation functions were calculated up to 1800 time steps but are only displayed to 1000 time steps so that the oscillatory character of the correlation functions shows up clearly. Also note the slower decay of these optical correlation&as compared to the INS correlation functions of fii. 1.

170

Fig. 3. (a) The OH-stretching region of the INS spectrum of fig. Id. The arrows indicate the experimentally measured position and width of the INS spectrum at 258 K from the data in ref. [3]. (b) The reduced infrared spectral intensity at 250 K calculated by Fourier transforming the IR correlation function of fig. 2a. The experimental infrared spectrum at 274 K [23] is shown for comparison. (c) The isotropic component of the reduced Raman spectral intensity at 250 K calculated by Fourier transforming the isotropic Raman correlation function of fig. 2b. The spectrum is plotted as 45Zi, in order to directly compare with the corresponding experimental data at 263 K [21]. (d) The anisotropic component of the reduced Raman spectral intensity at 250 K calculated by Fourier transforming the anisotropic Raman correlation function of fig. 2c. The spectrum is plotted as 41amso in order to compare directly with the corresponding experimental data at 263 K [21]. Note that the same intensity scale is used for (c) and (d).

Volume 132, number 2

CHEMICAL PHYSICS LETTERS

S December 1986

relation functions. The main reason for this difference between optical (IR and Raman) and INS correlation functions is that the IR and Raman correlation functions are calculated only for the OH stretching motion (i.e. the motion of the proton relative to that of the oxygen) and thus reflect the damped oscillator motion of the OH bond. On the other hand, the INS correlation function is calculated for all the motions of the proton and hence it has a different time dependence. The different behaviour of these time correlations is also due to the fact that the INS correlation function when expanded in powers of &r(t) is a sum of all orders in &r(t) whereas the Raman and IR time correlation functions involve only low-order terms when expanded in powers of &r(t). Fig. 3 gives the corresponding spectral functions, in the high-frequency OH-stretching region, as defined in eq. (8) for the reduced infrared and the isotropic and anisotropic Raman intensities. Also shown is the Q = 8.72 A-’ INS spectrum in the OH-stretch region. The experimental infrared spectrum [23], and the isotropic and anisotropic Raman spectra [21] are also shown on the same graph. These curves taken individually are qualitatively similar to their corresponding experimental band shapes. The peak positions in all three calculated spectra occur somewhat lower than observed experimentally. To some extent this is to be expected since the calculated spectrum corresponds to a deeply supercooled state and peak frequencies decrease as temperature decreases [24]. Also the calculated spectra are somewhat more structured than the experimental counterparts. This is also partly explained on the basis of supercooling and partly due to the fact that the simulation deals with a much smaller system and thus decreases the extent of configurational averaging. Also it should be recognized that the vibrational potential parameters were taken from the free molecule frequencies. For example, we have found that a small change in the parameter p from 2.37 to 2.566 A-’ shifts the peak from 3100 cm-’ to 3500 cm-‘. Clearly these parameters can be optimized to generate the correct frequencies for liquid water. The calculated band shapes for the infrared and neutron spectra are in good qualitative agreement with the experimentally measured spectra. The calculated Raman spectrum, in particular the isotropic component shows the largest deviation from the experiment. The calculated spectrum accentuates the intensity in the low-frequency side (below 3200 cm-l) relative to the intensity above 3200 cm-’ . This may partly be due to the fact that the spectrum shown in fig. 3c does not include Fermi resonance terms which contribute to the spectral intensity in this region. Furthermore there are undoubtedly second-order effects due to dipole-dipole coupling terms. The fact that the agreement is poorer for the isotropic component than the anisotropic component suggests that dipole-dipole coupling terms contribute more to the diagonal elements of the polarizability tensor than the off-diagonal terms. It should also be emphasized that the calculated spectrum was generated at 250 K and thus should be compared with the measured spectra in the supercooled state. The data of d’Arrigo et al. [24] do indeed show that the ratio of the spectral intensity of the low-frequency component relative to that of the high-frequency component increases with decreasing temperature. They also find that the peaks are shifted to lower frequencies as temperature decreases. We have confirmed these temperature-dependent effects by calculating the Raman spectrum at 325 K. These results will be discussed in detail elsewhere. It is interesting to note that Scherer et al. [21] used three peaks to fit the broad experimental spectrum. The three peaks that are well resolved in the calculated isotropic Raman spectrum at 3120,3200, and 3450 cm-’ correspond to the fitted peaks at 3225,342O and 3615 cm -’ in Scherer et al. [21]. It is interesting to note that in our model the infrared and the anisotropic Raman spectral band shapes are qualitatively very similar to the high4 inelastic neutron spectrum, whereas the isotropic Raman spectrum is quite different. Perhaps this similarity is related to some basic property of the density of states and should be looked into further. In conclusion, we note from a comparison of the three spectra that although the frequency range of the stretching bands is nearly identical in the different spectra, the band shapes are quite different from each other. The calculated spectra reproduce the broad band shape observed experimentally in liquid water quite well considering the fact that the potential parameters were in no way fitted to liquid state data. We therefore conclude that in “SPCflexible molecule” water the anharmonicity of the intramolecular potential and the coupling of the simulated intermolecular H-bond interactions are strong enough to give a stretching band shape which substantially deviates from the free molecule behaviour and do indeed show most of the liquid water spectral features.

171

CHEMICAL PHYSICS LETTERS

Volume 132. number 2

5 December 1986

Acknowledgement The reserach

of RB was supported

port from NSF. TB acknowledges

joint financial

by a grant from the Office of Naval Research,

a grant-in-aid

of research

SHC acknowledges

partial

from BU chapter of Sigma Xi. MAR acknowledges

support of NATO and CNR.

References [l] [2] [3 ] [4] (51 [6] [7] [8] [9] [lo] [ll]

J.R. Scherer, in: Advances in infrared and Raman spectroscopy, Vol. 5 (Heyden, London, 1978) ch. 3. G.E. Walrafen, in: Water: a comprehensive treatise, Vol. 1, ed. F. Franks (Plenum Press, New York, 1971) ch. 5. S.H. Chen, K. Toukan, C.K. Loong, D.L. Price and J. Teixeira, Phys. Rev. Letters 53 (1984) 1360. S.A. Rice, M.S. Bergren, A.C. Belch and G. Nielson, J. Phys. Chem. 87 (1983) 4295. A.C. Belch and S.A. Rice, J. Chem. Phys. 78 (1983) 4817. J.R. Reimers and R.O. Watts, Chem. Phys. 91 (1984) 201, and references therein. J.P.M. Postma, H.J.C. Berendson and T.P. Straatsma, J. Phys. (Paris) 45 (1984) C7. H.L. Lemberg and F.H. Stillmger, J. Chem. Phys. 62 (1975) 1677. F.H. Stillinger and C.W. David, J. Chem. Phys. 73 (1980) 3333, and references therein. K. Toukan and A. Rahman, Phys. Rev. B31 (1985) 2643. H.J.C. Berendson, J.P.H. Postma, W.F. van Gunsteren and J. Hermans, in: Intermolecular forces, ed. B. Pullman (Reidel, Dordrecht, 1981). [12] D.R. Fredkin, A. Komornicki, S.R. White and K.R. Wilson, J. Chem. Phys. 78 (1983) 7077. [13] A. Rahman and F.H. Stlllinger, Phys. Rev. A10 (1974) 368. [14] A.K. Soper, Chem. Phys. 88 (1984) 187. [15] S.H. Chen and S. Yip, Physics Today,29 (1976) 32. [16] M.A. Ricci, S.H. Chen, D.L. Price, C.K. Loong, K. Toukan and J. Teixeira, Physica B136 (1986) 190. [ 171 A.C. Lie and E. Clementi, Phys. Rev. A33 (1986) 2679. [18] D.A. McQuarrie, Statistical mechanics (Harper and Row, London, 1976). [19] S. Ikawa and S. Maeda, Spectrochim. Acta A24 (1968) 655. [20] J.R. Scherer and R.G. Snyder, J. Chem. Phys. 67 (1977) 4794. [21] J.R. Scherer, M.K. Go and S. Kint, J. Phys. Chem. 78 (1974) 1304. 1221 B.J. Berne and R. Pecora, Dynamic light scattering (Wiley, New York, 1975). [23] L.W. Pinkley, P.P. Sethna and D. Williams, J. Opt. Sot. Am. 67 (1977) 494. [24] G. d’Arrigo, G. Maisano, F. Mallamace, P. Migliardo and F. Wanderhngh, J. Chem Phys. 75 (1981) 4264.

172

sup-

the