Role of dipolar correlations in the IR spectra of water and ice
We report simulated infrared spectra of deuterated water and ice using Car-Parrinello molecular dynamics with maximally localized Wannier functions. Experimental features are accurately reproduced within the harmonic approximation. By decomposing the lineshapes in terms of intra and intermolecular dipole correlation functions we find that short-range intermolecular dynamic charge fluctuations associated to hydrogen bonds are prominent over the entire spectral range. Our analysis reveals the origin of several spectral features and identifies network bending modes in the far IR range.
pacs:36.20.Ng, 71.15.Pd, 32.10.Dk, 78.20.Ci
It is well known that correlations among the dipole moments of hydrogen bonded molecules greatly enhance the static dielectric response of water and ice (see e.g. Ref. Sharma07 () and references therein). However, the role played by dipolar correlations in the dynamic dielectric response of hydrogen bonded systems is not well understood, although it has often been speculated that intermolecular couplings affect significantly infrared (IR) and Raman spectra Whalley77 (); Rice83 (); Buch99 (); Guillot91 (); Silvestrelli97 (). A difficulty in modeling these spectra is that the dynamic dielectric susceptibility depends on the adiabatic response of the electrons to nuclear dynamics. Another difficulty, in principle even greater, is that nuclear quantum effects cannot be discarded a priori in water systems at, or below, room temperature, and it is beyond current theoretical and computational capabilities to fully account for dynamic quantum effects in systems of this complexity. Recent studies have shown, however, that the IR absorption lineshapes of liquid water and ice are overall reproduced remarkably well by ab initio molecular dynamics (AIMD) simulations, which treat nuclear dynamics classically but include explicitly, and accurately, the quantum adiabatic response of the electrons Sharma05 (); Iftimie05 ().
In this paper we adopt the AIMD framework and use maximally localized Wannier functions (MLWF) to separate intra and inter molecular contributions to the dynamic dipolar correlations in liquid water and ice. We find that hydrogen bonds (H-bonds) make short-range intermolecular fluctuations of the electronic charge as important as the intramolecular fluctuations over the entire spectral range. Our analysis sheds new light into the origin of several spectral features and provides an unambiguous identification of network bending modes in the far IR region Zelsmann95 ().
The IR absorption coefficient per unit length is related to the refractive index and the imaginary part of the dielectric constant by . Within linear response theory, is given by the power spectrum of the time correlation function of the total dipole operator McQuarrie00 (). Following common practice, we approximate the quantum time correlation function with the classical one, i.e. with , where is the total dipole moment in the simulation cell and the brackets indicate classical ensemble average. Since there is only one classical correlation function while the quantum time correlation function can be expressed in several equivalent ways, this substitution leads to formulae for the IR absorption coefficient characterized by different prefactors known as quantum correction factors Ramirez04 (). We adopt here the so-called harmonic approximation (HA) which follows by replacing the Kubo-transformed quantum correlation function with the classical one, leading to:
where is the inverse temperature. In the harmonic regime HA is exact Bader94 (). Ramírez et al. Ramirez04 () showed that HA is the only correction factor that satisfies the fluctuation-dissipation theorem in addition to detailed balance. The same authors also found that HA performs better than the other quantum correction factors for one dimensional anharmonic potentials that model different H-bond scenarios.
Ii Computational details
In this work, we compute the IR spectra of DO ice and water using a Car-Parrinello (CP) scheme Car85 (), in which maximally localized Wannier functions (MLWFs) Marzari97 () represent “on the fly” the electronic wave functions Sharma03 (). In the simulation of water (ice), we use a periodically repeated cubic (orthorhombic) cell with 64 (96) DO molecules at the density of 1.1 (1.0) g/ml. For ice we use a proton disordered ice Ih configuration generated as in Ref. Sharma07 (). We adopt norm-conserving pseudopotentials Troullier91 () and the DFT Perdew-Burke-Ernzerhof (PBE) functional for exchange and correlation Perdew96 (), using a plane-wave cutoff of 60 Ry for water and of 85 Ry for ice Details_Cutoff (). We use an integration time step of 7 a.u. (0.17 fs) and a fictitious electron mass of 350 a.u. in the CP equations Car85 (); Sharma05 (); Sharma07 (). The temperature is controlled by coupling the nuclei to a single Nosé thermostat Blochl92 () with mass . In this way the average temperature of the system is set to 330 K in water and to 268 K in ice.
At each time step we assign a dipole to the -th molecule in the cell. Here are the positions of the nuclei ( = D, D, O) and of the four MLWF centers (, with ) associated to the eight valence electrons of a water molecule Silvestrelli99 (); Sharma05 (); Sharma07 (). The -th dipole is conventionally located at the molecular center of mass , and depends on the local environment. The total dipole moment is Pasquarello03 (); Resta94 () and the IR spectra are computed with Eq. (1). The time correlation function is averaged over an equilibrated MD trajectory of 10 ps in ice and of 23 ps in water. A Gaussian window function Harris1978 () is used in the discrete Fourier transform.
Iii Results and discussion
iii.1 Calculated IR spectra
As in previous HA calculations Sharma05 (); Iftimie05 () the overall agreement between theory and experiment is good. For ice, the positions of the band maxima corresponding to H-bond stretching, libration, bending, combination, and OD stretching are: 200 (222), 660 (640), 1140 (1210), 1660 (1650), and 2120 (2425) cm respectively (wavenumbers in parentheses are experimental values from Ref. Eisenberg69 ()). For water, the corresponding values are: 200 (187), 510 (505), 1140 (1215), 1550 (1555), and 2160 (2450) cm respectively. The details of the experimental features are well reproduced, e.g. the asymmetric shape of the OD stretching bands, the skewed ice libration band, and even the small combination bands. The good performance of the harmonic approximation may appear somewhat surprising given the large shifts and broadenings of the vibrational frequencies of an isolated molecule that arise in condensed phase. Furthermore, significant anharmonicity has been detected in the excited hydrogen stretching modes in liquid water Bakker2002 (). IR absorption, however, probes equilibrium properties which are dominated by the vibrational ground-state in the hydrogen stretching region. Then, the overall similarity between calculated and observed spectra suggests that the main effects of anharmonicity are sufficiently well captured by classical dynamics.
The largest discrepancy between simulation and experiment occurs in the hydrogen stretching modes, which are redshifted compared to experiment by 300 cm. Part of this error originates from our choice of the mass for the fictitious dynamics of the electrons. We can quantify this effect by calculating the IR stretching band in ice from Born-Oppenheimer (BO) molecular dynamics simulations, in which the electrons are kept in the instantaneous ground-state by minimizing the energy functional at each time step. In CP simulations this condition is enforced dynamically and the outcome depends on the choice of . BO calculations show that a redshift of up to 80 cm in the IR stretching band of ice can be attributed to our choice of Details_mu (). However, even when the BO separation is strictly enforced the IR stretching band in ice is more than 200 cm below experiment. The adopted DFT approximation is a likely cause of this error. For instance, in the HO water monomer, PBE yields stretching frequencies that are 140 cm lower than the corresponding experimental values Xu2004 (). Taking the isotopic mass effect into account we should expect a redshift of 100 cm for the stretching band of DO in gas phase. The actual effect that we observe in condensed phase, 200 cm, is larger than that and is consistent with the known H-bond over binding present in PBE water. Other sources of inaccuracy in the calculated spectra include the finite basis set, the cell size and the effect of temperature. Quantum anharmonic corrections are beyond our approach.
The similarity between the spectra of solid and liquid water is consistent with the standard picture of a tetrahedral, ice-like local structure of the liquid. With more disorder in the liquid phase, modes close in frequency overlap with each other, leading to broader bands than in the solid. Fluctuating local fields and coupling between neighboring molecules cause additional broadening. The relative importance of local fields and intermolecular couplings can be quantified in terms of intra and inter spectral contributions.
iii.2 Intra vs. inter contribution
While separating intra and intermolecular contributions would be difficult or even impossible in experiment, in theory it can be easily achieved due to the additive nature of two-body correlation functions. In Fig. 2 we report the intra and inter contributions to the IR spectra of ice and water.
The total spectrum is computed from Eq. (1), in which . The intra contribution corresponds to retaining only the terms with , the remaining terms add up to the inter contribution.
In water, the most prominent band around 2100 cm is due to the overlap of the symmetric and asymmetric OD stretching modes together with the overtone of the bending mode Eisenberg69 (); Whalley77 (). With two maxima, the inter contribution has larger intensity and lower frequency compared to its intra counterpart. This gives rise to the three broad features and the asymmetric shape of the band observed in experiments. Similar effects occur in ice, where three main features have been experimentally identified in the stretching band Bergren78 (); Buch99 (). In this case, however, we can not identify distinct intermolecular features within our chosen Fourier filter Details_smooth (). Two peaks resembling the experimental features appear in the inter contribution when we reduce the filtering but they are too closely separated compared to experiment. In contrast to the stretching modes, the inter contributions to the bending modes around 1140 cm are negative in both ice and water. This effect has a simple intuitive explanation: due to the electrostatic attraction between opposite charges, OD stretching modes are correlated between neighboring molecules whereas bending modes are anti-correlated. In addition, we find that the broad shape of the small combination band at 1550 cm in water (1660 cm in ice) is mostly due to the inter contribution.
The far infrared region is more complex. Two features have been identified in the libration band of the Raman spectrum of the liquid, but only one mode was believed to be IR active on the basis of a simple tetrahedral model with C symmetry Eisenberg69 (). A more recent IR experiment, however, suggested that both modes are visible in the liquid Zelsmann95 (). These modes have their counterpart in ice and, indeed, we find two IR features at 460 cm and at 660 cm in our calculated spectrum. Both features have intra and inter contributions but, interestingly, the inter contribution at 460 cm is negative, largely cancelling its intra counterpart. Instead the inter contribution at 660 cm is small and positive, slightly enhancing the corresponding spectral feature. As a result, the ice libration band has the skewed shape observed in experiments (Fig. 1). In the liquid, the situation is similar with the features shifted to 370 cm and 590 cm. Broadening due to disorder makes the lineshape look more symmetric.
The features at 50 cm in water and at 200 cm in water and ice are H-bond network modes with bending and stretching character, respectively Zelsmann95 (); Sharma05 (). In Ref. Sharma05 (), the inter character of the 200 cm mode was emphasized. A more accurate analysis reveals, however, that both intra and inter contributions are present, with an inter/intra ratio of 2:1 in ice and of 1:1 in water. This analysis is facilitated in the case of ice, where the network modes are well separated from the libration band.
iii.3 Identification of the network bending mode
Interestingly, we also understand the origin of the feature at 50 cm. Albeit weak, this feature has been identified in experiments particularly below room temperature Zelsmann95 (). In our simulation it appears as a weak intramolecular feature in both water and ice, where the peak is at 70 cm. Its intra character is surprising because the feature originates from network bending modes, as shown by the power spectrum of the autocorrelation function of the angle between two adjacent H-bonds, i.e. the angle defined by the oxygen atoms of three molecules linked by H-bonds. The corresponding spectrum for ice is reported in the inset of Fig. 2 and is sharply peaked at 70 cm. We also report in the same inset the power spectrum of the autocorrelation function of the angle between the two covalent OD bonds in a molecule. As expected this spectrum is sharply peaked at 1140 cm but, interestingly, it also exhibits a weaker feature at 70 cm due to the modulation of the molecular angle induced by the network bending modes. Thus, the absence of inter character in the IR spectrum derives from the negligible intermolecular charge fluctuation associated to network bending modes, while the intra contribution reflects the modulation of the molecular dipole moment induced by these modes.
iii.4 Spatial extent of dynamic dipolar correlations
To gain insight on the spatial extent of the intermolecular correlations, we extend the approach of Ref. Sharma07 () to the dynamic domain. We define the density of the intermolecular dipole correlation function as:
where is the distance between the center of mass of molecule at time and molecule at time . if at time molecule is in spherical shell centered on molecule at time , whereas otherwise. The corresponding integrated contribution to the IR spectrum up to distance is given by:
The dependence of several peak intensities is reported in Fig. 3.
The curves indicate that the most important intermolecular correlations for all the modes occur when is comprised between 4 and 6 bohr, i.e. in correspondence with the first coordination shell, in both water and ice. Beyond 6 bohr only minor additional contributions are present, similar to the findings for static correlations in Ref. Sharma07 (). Short-range dynamic correlations are a direct manifestation of the presence of H-bonds between adjacent molecules. While static dipolar correlations due to H-bonds always enhance the dielectric response, dynamic correlations do enhance the system response at some frequencies but they also suppress it at other frequencies.
Given the importance of intermolecular correlations, only models that include their effect at a sufficient level of detail can realistically describe the observed experimental features. An important result of our study is that these correlations play an equally important role over the entire spectral range and are not limited to the network modes that characterize the far IR range. We expect that similar dynamic correlations should also be crucial to interpret other spectroscopic data, such as Raman, which probe the adiabatic dynamics of the electrons. It would be difficult to include properly these effects in simulations based on classical empirical potentials. Our analysis shows that liquid water and ice are systems in which the molecules are strongly correlated, suggesting that their IR spectra cannot be simply explained in terms of a weakly perturbed collection of individual molecules. Finally, the decomposition technique introduced in this paper is not restricted to bulk water and ice, but can also be applied to other systems and to water at interfaces. In the latter case one may expect that characteristic signatures of hydrophilic and hydrophobic interfaces should emerge from this analysis.
Acknowledgements.We would like to thank Davide Donadio for useful discussions. W.C. and R.C. gratefully acknowledge support from the NSF-MRSEC grant DMR-0213706. M.S. and G.G. gratefully acknowledge support from Scidac grant No. DE-FG02-06ER46262.
- (1) M. Sharma, R. Resta, and R. Car, Phys. Rev. Lett. 98, 247401 (2007).
- (2) E. Whalley, Can. J. Chem. 55, 3429 (1977).
- (3) S. A. Rice, M. S. Bergren, A. C. Belch, and G. Nielsen, J. Phys. Chem. 87, 4295 (1983).
- (4) V. Buch and J. P. Devlin, J. Chem. Phys. 110, 3437 (1999).
- (5) B. Guillot, J. Chem. Phys. 95, 1543 (1991).
- (6) P. L. Silvestrelli, M. Bernasconi, and M. Parrinello, Chem. Phys. Lett. 277, 478 (1997).
- (7) M. Sharma, R. Resta, and R. Car, Phys. Rev. Lett. 95, 187401 (2005).
- (8) R. Iftimie and M. E. Tuckerman, J. Chem. Phys. 122, 214508 (2005).
- (9) H. R. Zelsmann, J. Mol. Struct. 350, 95 (1995).
- (10) D. A. McQuarrie, Statistical Mechanics (University Science Books, Sausalito, California, 2000).
- (11) R. Ramírez, T. López-Ciudad, P. K. P., and D. Marx, J. Chem. Phys. 121, 3973 (2004).
- (12) J. S. Bader and B. J. Berne, J. Chem. Phys. 100, 8359 (1994).
- (13) R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
- (14) N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 (1997).
- (15) M. Sharma, Y. Wu, and R. Car, Int. J. Quantum Chem. 95, 821 (2003).
- (16) N. Troullier and J. L. Martins, Phys. Rev. B 43, 1993 (1991).
- (17) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- (18) Increasing the plane wave cutoff in the water simulation to 85 Ry would not change our qualitative results. The main effect would be variations of a few percent in the peak frequencies.
- (19) P. E. Blöchl and M. Parrinello, Phys. Rev. B 45, 9413 (1992).
- (20) P. L. Silvestrelli and M. Parrinello, Phys. Rev. Lett. 82, 3308 (1999).
- (21) A. Pasquarello and R. Resta, Phys. Rev. B 68, 174302 (2003).
- (22) R. Resta, Rev. Mod. Phys. 66, 899 (1994).
- (23) F. J. Harris, Proc. IEEE 66, 51 (1978).
- (24) In the calculated IR spectra in Fig. 1 and 2 a Fourier filter is used to eliminate fluctuations with a width smaller than 60 cm.
- (25) J. E. Bertie, H. J. Labbé, and E. Whalley, J. Chem. Phys. 50, 4501 (1969).
- (26) M. S. Bergren, D. Schuh, M. G. Sceats, and S. A. Rice, J. Chem. Phys. 69, 3477 (1978).
- (27) J. E. Bertie, M. K. Ahmed, and H. H. Eysel, J. Phys. Chem. 93, 2210 (1989).
- (28) D. S. Eisenberg and W. Kauzmann, The Structure and Properties of Water (Clarendon, Oxford, 1969).
- (29) H. J. Bakker and H.-K. Nienhuys, Science 297, 587 (2002).
- (30) The shift in the molecular frequencies induced by the fictitious dynamics of the electrons can be reduced by lowering the value of at the price of increasing the computational cost. For instance, using a.u. would double the computational cost. Calculations on the isolated water monomer show that the CP shift in the stretching frequencies is essentially negligible when is smaller than 100 a.u.. Notice that the influence of the fictitious electron dynamics on the nuclear dynamics decreases rapidly with the frequency separation between electronic and nuclear motions. As a consequence the CP redshift is larger on the nuclear stretching modes which have frequency closer to the minimum electronic frequency. To very good approximation the latter is given by in terms of the Kohn-Sham band gap Car2006 ().
- (31) R. Car, in Conceptual Foundations of Materials: A Standard Model for Ground- and Excited-State Properties, edited by S. G. Louie and M. L. Cohen (Series: Contemporary Concepts of Condensed Matter Science, Elsevier, Amsterdam, 2006), Chap. 3, p. 64.
- (32) X. Xu and W. A. Goddard, III, J. Phys. Chem. A 108, 2305 (2004).