J. Phys. Chem. 1996, 100, 1411-1414
1411
Frequency- and Wave-Vector-Dependent Susceptibility of Water† Andrij Trokhymchuk*,‡ Instituto de Quimica de la UNAM, Circuito Exterior, Coyoacan 04510, Mexico DF
Myroslav Holovko Institute for Condensed Matter Physics, National Academy of Sciences of the Ukraine, LViV 11, Ukraine ReceiVed: August 15, 1995; In Final Form: October 20, 1995X
The charge susceptibility function of a flexible water model is investigated using the theory of the dynamical dielectric function of interaction-site model solvent developed by H. L. Friedman and his colaborators. The atom-atom pair correlation functions used as input in this approach were derived from results of molecular dynamic computer simulations. They were extended to larger separations by the use of integral equation theory techniques. It is pointed out that a more careful description of long-ranged correlations can provide new information about the nature of the collective dynamics in water.
Investigation of the dynamical properties of water is one of the more challenging areas in the modern statistical-mechanical theory of solutions,1-9 as well as experimental studies using ultrafast spectroscopic techniques.10,11 A major problem is the lack of agreement between the theoretical and experimental results; therefore, several interesting key questions still remain unresolved.12 The reason for this is not fully clear. From the theoretical point of view there are two main sources of error: a more accurate theory is needed, and a more realistic model assumption should be used. The widely used approximation in the analytical theories reported until now is the assumption on a linear response of a solvent to change in the charge distribution of a solute. The simplest model of a solvent studied extensively is based on a dielectric continuum description.2,3 The most recent developments have focused on the solvent as an assembly of the dipolar hard sphere molecules.4,9 In ref 1 the molecular theory of the frequency-dependent and wave-vector-dependent dielectric properties of an interactionsite model solvent has been derived. In the framework of this theory the dynamical properties of a polar solvent are expressed in terms of the charge susceptibility function χµφ0(k,ω) which is associated with the dynamics of a inhomogeneous solvent in the presence of a solute in its initial charge state. The molecular expression for the external charge susceptibility has the form1
[
Cµµ(k,ω)
χµφ0(k,ω) ) βSµµ(k) iω
Sµµ(k)
-1
]
(1)
where Sµµ(k) is the charge-charge static polarization structure factor which for an interaction-site model may be expressed in terms of the atom-atom intra- and intermolecular correlation functions.13 Cµµ(k,ω) is the Fourier-Laplace transform of the time correlation function Cµµ(k,t) of the spatial Fourier component of the microscopic charge density
Fˆµ(k) ) ∑[∑qjeik‚δraj]eik‚ra a
(2)
j∈a
Here ra is the position of the center of mass of molecule a, δraj is the position of atom j in molecule a, and qj is its partial charge. †
Dedicated to Harold L. Friedman on the occasion of his 70th birthday. Permanent address: Institut for Condensed Matter Physics, National Academy of Sciences of the Ukraine, Lviv 11, Ukraine. X Abstract published in AdVance ACS Abstracts, December 1, 1995. ‡
0022-3654/96/20100-1411$12.00/0
The reference memory function approximation was applied in ref 14 to formulate a useful approximation to C ˆ µµ(k,t), which gives
Cµµ(k,ω)
iω
Sµµ(k)
)
N(ω) - 1 N(ω)[1 - 0R(k)] - 1
(3)
where the function N(ω) ) ((ω) - 1)/(0 - 1) is related to the dielectric permittivity 0 and the frequency-dependent dielectric function (ω) and, particularly, can be determined by computer experiment for the corresponding water model;15 R(k) ) RJ(k)/ RS(k) is the renormalization static factor defined as a ratio of the kinetic part, RJ(k) ) Jµ,L(k)/Jµ,L(k)0), and the structural part, RS(k) ) {Sµµ(k)/k2}/{Sµµ(k)/k2}k)0, where Jµ,L(k) is the value of the time correlation function Jµ,L(k,t)0) of the longitudinal part of the charge current density Fˆ˙ µ(k)/(ik). Thus, the above scheme allows calculation of the dynamical properties for interaction-site models of water; the frequency-dependent dielectric function (ω) and the static charge-charge polarization structure factor Sµµ(k) serve as input in this theory. Recently, it has been applied14 to calculate the dynamical properties of TIP4P and SPC/E water using Neumann’s parametrization15 of (ω), and the extended reference interactionsite method (XRISM), to calculate Sµµ(k). It was shown that, due to the limited capability of XRISM to reproduce the structure of model water, the results are sensitive to a model assumption that influences mainly the renormalization factor R(k). XRISM theory reproduces the value {Sµµ(k)/k2}k)0 correctly but only approximates the behavior of this function in the intermediate region of wave vectors k, which is important in the understanding of the collective dynamics of water. In refs 13, 16, and 17 we have used an approach18 based on the combination of computer simulation and integral equation theory techniques to calculate the small wave numbers behavior of a static water structure factor. The computer simulation result employed is an MD simulation run18 with the BJH flexible version of the central-force potential functions19 which was performed for 200 water molecules at 292 K and a density of 0.9718 g/cm3 giving a side length of the basic periodic cube of 18.17 Å. The simulated radial distribution functions have been extended to 25.6 Å, as is described in ref 18. The distance 25.6 Å follows from the procedure of numerical solution of the Ornsten-Zernike equations with 1024 grid points at a spacing © 1996 American Chemical Society
1412 J. Phys. Chem., Vol. 100, No. 4, 1996
Trokhymchuk and Holovko
of 0.025 Å. The radial distribution functions gij(r) calculated in such a way are discussed in refs 17 and 18. Particularly, the first two minima and the second maximum calculated for gOO(r) are slightly more pronounced, though the first peak agrees well with the scattering data. At large distances (r > 12 Å), all three gij(r) oscillate around unity with a period significantly larger than that at small and intermediate distances. Starting from the distances r ) 13.3, 13.2, and 12.7 Å, the oxygenoxygen, hydrogen-hydrogen, and oxygen-hydrogen radial distribution functions, respectively, are less than unity up to r ≈ 20 Å, where they become larger than unity. At larger distances they approach unity monotonically again. The charge-charge polarization structure factor Sµµ(k) for this model can be represented as13
4πβSµµ(k) ) Sµµ(k)intra + Sµµ(k)inter
(4)
where intra (k) ) Sµµ
9y 3 1 + wHH(k) - 2wOH(k) L2 2 2
[
]
(5)
9y ∫[gOO(r) + gHH(r) - 2gOH(r)]sinkrkrr2 dr L2 (6)
inter Sµµ (k) ) 4πF
with
wij(k) )
sin kLij -Bij2k2 e kLij
(7)
the Fourier transform of the intramolecular atom-atom correlation functions, and y ) (4/9)πFβµ2, L2 ) LOH2 - 0.25LHH2. In the last expressions, µ is a permanent dipole moment, and Lij, Bij are the average values of the bond lengths and the associated root-mean-square deviations of the water molecules, respectively. For the BJH model these have been calculated:13 LOH ) 0.975 Å, LHH ) 1.500 Å, BOH ) 0.032 Å, and BHH ) 0.063 Å. In particular, we found that Sµµ(k) has an oscillating behavior in the region of small k which yields18 well-resolved minima and maxima of the calculated diffraction intensity around 0.31 and 0.62 Å-1, respectively. However, the oscillating behavior of measured scattering intensity at small k, in fact, has not been discussed in diffraction studies of water (the presence of minima around 0.4 Å-1 has been observed20,21). We suppose this is not an artifact introduced by the BJH model since this water model has proved its usefulness in various simulations of water and electrolyte solutions.22 Nevertheless, it cannot be excluded that the hydrogen bond strength in the BJH model might be exaggerated, and the magnitudes of the minima and maxima mentioned are in reality less pronounced than those following from model calculations. We have also shown16 that the behavior of the structure factor at small and intermediate values of wave vector is bound up with that of radial distribution functions at separations larger than 10 Å. Obviously, such behavior of the structure factor does not come from computer simulation results obtained via Fourier transform of radial distribution functions. At the same time, calculating the structure factor of ST2 water directly from computer simulation23 gives the diffraction intensity behavior similar to that of BJH water.17 In ref 16 we discussed the relation between such behavior of structure factors and the existence of hydrogenbonded formations in water. Applying Sµµ(k) to calculate the static dielectric properties13 leads to an unusual splitting of the negative branch of the static dielectric function L(k) into two parts in this region. It indicates the existence of different types
Figure 1. Renormalization static factor R(k) for the BJH water model. The long-dashed line is the structural part RS(k), the short-dashed line is the kinetic part RJ(k), and the solid line is R(k).
of contributions to the polarization fluctuations: one of them is related to correlation between monomeric molecules, while a second one is due to correlation between those involved in hydrogen bonds. It causes the different screening properties of water: from “perfect” screening (negative values of L(k)) up to “complete” screening (large positive values of L(k)). Simple calculations17 of the dispersion curve ω(k) in the framework of the undamped theory within the mentioned reference memory function approximation qualitatively describe the experimental data obtained from inelastic neutron scattering measurements,24 which demonstrates the presence of damped density fluctuations propagating with a velocity larger than twice the hydrodynamic sound velocity. The present work addresses application of Sµµ(k) in the form of eqs 4-7 to investigate the frequency-dependent and wavevector-dependent susceptibility function of a flexible water model at room temperature. Figure 1 shows the quantities RJ(k), RS(k), and R(k) calculated for the BJH model. Since Sµµ(k) is calculated17 in the region k e 0.62 Å-1 with an increment of 0.31 Å-1, we calculate the structural part RS(k) in the same way. The method for calculating the kinetic part RJ(k) for the BJH model is similar to that for the SPC/E model described in the Appendix of ref 14. A new feature of R(k) caused by the oscillating behavior of Sµµ(k) at low values of k is the presence of a narrow region k where R(k) changes its sign. This region coincides with that of the positive branch in the longitudinal static dielectric function L(k) discussed in refs 13 and 17. Since the frequency-dependent dielectric function (ω) was not determined for the BJH model, we use here Neumann’s parametrization15 for TIP4P water at 293 K with some changes. Namely, we replace the value of dielectric permittivity 0 ) 52.6 in this formula by that for the BJH model, 0 ) 79.29.13 We also replaced parameters ω0 f (ω02 + γ2)1/2 and tan Φ f -γ/(ω02 + γ2)1/2. The frequency-dependent dielectric constant (ω) ) ′(ω) - i′′(ω) calculated in this way is displayed in Figure 2. As shown with comparision with experimental data the model (ω) satisfactorily describes a wide frequency region up to approximately 400 ps-1. The real and imaginary parts of the charge susceptibility χµφ0(k,ω) calculated according to (1) are depicted in Figure 3 at some selected values of wave vector k ) 0, 0.31, 0.62, 0.8, 1.5, and 3.0 Å-1. At k ) 0 there are two well-resolved peaks: in the high-frequency region very close to the frequency of librational resonance in (ω) and in the low-frequency region associated with diffusional motions. Similar behavior was
Susceptibility Function of Flexible Water Model
J. Phys. Chem., Vol. 100, No. 4, 1996 1413
Figure 2. Real and imaginary parts of the frequency-dependent dielectric constant (ω) of water at room temperature. Solid lines are the calculation for the BJH model; circles are the experimental results taken from ref 15. Frequency is shown in ps-1. The imaginary part has been shifted by a factor of 0.1.
Figure 4. Function ω2Sµµ(k,ω)/k2 of BJH water. Each subgraph has the same ordinate scale from 0 to 70.
To relate the behavior of the susceptibility function with the observed features of water, in Figure 4 we display the behavior of the function ω2Sµµ(k,ω)/k2 for the same values of k. The charge-charge dynamic structure factor Sµµ(k,ω) is calculated by14
Sµµ(k,ω) )
Figure 3. Real (solid line) and imaginary (dashed line) parts of the susceptibility function χµφ0(k,ω). Frequency is shown in ps-1.
observed26 for the dynamic dielectric response function of water at k f 0 constructed by combining the dynamical self-part of the van Hove correlation function with a static structural information from the hypernetted chain integral equation theory. These two peaks persist with an increase in k. In the region where we found an oscillating behavior of the diffraction intensity, the diffusion peak significantly increases and is slightly shifted to lower frequencies; the high-frequency peak remains virtually unchanged. In the narrow region of k values around 0.8 Å-1 where the longitudinal static dielectric function L(k) is positive17 (indicating increase in the screening), the behavior of the susceptibility function is qualitatively different (indicating decrease in the response): diffusion peaks are absent, and only very small peaks shifted to high frequencies are observed. As k increases (k approximately 1.1 Å-1), the response increases, and the susceptibility function displays the behavior similar to that at k f 0. In the region of wave-vector values around 3 Å-1 where the static structure factor of water obtained from both theoretical calculations and measurements has a minimum bound up with splitting of the main peak, the susceptibility function has a form similar to that around 0.3-0.6 Å-1 with a high-diffusion peak.
1 Re{Cµµ(k,ω)} π
(8)
The observed peaks in Sµµ(k,ω) are situated in the same frequency region as the optical-like collective mode reported from molecular dynamic calculations27 and can be identified with high-frequency collective excitations in water. The observation of coherent inelastic neutron scattering for shortwavelength high-frequency collective excitations in heavy water24 indicates dramatic increases in damping at k > 0.5 Å-1, which we relate to the behavior of Sµµ(k,ω) at k around 0.8 Å-1. The partial structure functions derived from a combination of computer simulation and integral equation theory techniques are applied to investigate the dynamical susceptibility function in liquid water. In the limit k f 0 the behavior of the susceptibility function is similar to that reported in literature. At short and intermediate values of wave vector the susceptibility has new features caused by an oscillating behavior of the static structure factor in this wave-vector region. The behavior of the charge-charge dynamic structure factor calculated in the same approach shows the existence of high-frequency collective excitations with qualitatively the same wave-vector dependence as those observed by coherent inelastic neutron scattering. It follows that a more careful calculation of the radical distribution functions at large separations (or partial structure factors at small k) can provide new information about the nature of the collective dynamics of water molecules in a liquid state. Similar calculations for other water models, investigation of the temperature dependence of the result, as well as experimental studies in this field could clarify this problem. Acknowledgment. A.T. expresses his gratitude to the CONACYT of Mexico for the financial support of the project. M.H. acknowledges the support from the ISF Grant No. U1 J000. We are indebted to Dr. D. Williamson for critically
1414 J. Phys. Chem., Vol. 100, No. 4, 1996 reading the manuscript. Fruitful discussions with K. Heinzinger are acknowledged as well. References and Notes (1) Raineri, F.; Zhou, Y.; Friedman, H. L.; Stell, G. Chem. Phys. 1991, 152, 201. (2) Bagchi, B.; Castner, E. W.; Fleming, G. R. J. Mol. Struct. 1989, 194, 171. (3) van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 1985, 89, 4181. (4) Chandra, A.; Bagchi, B. J. Phys. Chem. 1990, 94, 1874. (5) Friedrich, V.; Kivelson, D. J. Chem. Phys. 1987, 86, 6425. (6) Calef, D. F.; Wolynes, P. G. J. Chem. Phys. 1983, 78, 4145. (7) Nichols, A. L., III; Calef, D. F. J. Chem. Phys. 1988, 89, 3783. (8) Loring, R. F.; Mukamel, S. J. Chem. Phys. 1987, 87, 1272. (9) Schwartz, B. J.; Rossky, P. J. J. Phys. Chem. 1995, 99, 2953. (10) Rosenthal, S. J.; Xie, X.; Du, M.; Fleming, G. R. J. Chem. Phys. 1991, 95, 4715. (11) Jimenez, R.; Fleming, G. R.; Kumar, P. V.; Maroncelli, M. Nature 1994, 369, 471. (12) Nandi, N.; Roy, S.; Bagchi, B. J. Chem. Phys. 1995, 102, 1390. (13) Trokhymchuk, A.; Holovko, M. F.; Heinzinger, K. J. Chem. Phys. 1993, 99, 2964. (14) Resat, H.; Raineri, F. O.; Friedman, H. L. J. Chem. Phys. 1992, 97, 2618. (15) Neumann, M. J. Chem. Phys. 1986, 85, 1567.
Trokhymchuk and Holovko (16) Trokhymchuk, A.; Holovko, M. F.; Heinzinger, K. Z. Naturforsch. 1995, 50a, 18. (17) Trokhymchuk, A.; Holovko, M. F.; Heinzinger, K. Mol. Phys., in press. (18) Trokhymchuk, A.; Holovko, M. F.; Spohr, E.; Heinzinger, K. Mol. Phys. 1992, 77, 903. (19) Bopp, P.; Jancso, G.; Heinzinger, K. Chem. Phys. Lett. 1983, 98, 129. (20) Bosio, L.; Teixeira, J.; Bellisent-Funel, M. C. Phys. ReV. A 1989, 39, 6612. (21) Xie, Y.; Ludwig, K. F., Jr.; Morales, G.; Hare, D. E.; Sorensen, C. M. Phys. ReV. Lett. 1993, 71, 2050. (22) Heinzinger, K. In Computer Modelling of Fluids, Polymers and Solids; Catlow, C. R. A., et al., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1990; p 357. (23) Stillinger, F. H.; Rahman, A. J. Chem. Phys. 1974, 60, 1545. (24) Teixeira, J.; Bellisent-Funel, M. C.; Chen, S. H.; Dorner, B. Phys. ReV. Lett. 1985, 54, 2681. (25) Streett, W. B.; Tildesley, D. J.; Saville, G. ACS Symp. Ser. 1978, 86, 144. (26) Kim, S.-H.; Vignale, G.; DeFacio, B. Phys. ReV. E 1994, 50, 4618. (27) Ricci, M. A.; Rocca, D.; Ruocco, G.; Vallauri, R. Phys. ReV. A 1989, 40, 7226.
JP952385R
