New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Hydrogen bonding and Raman, IR, and 2DIR spectroscopy of dilute HOD in liquid D_{2}O

Edited by Robin M. Hochstrasser, University of Pennsylvania, Philadelphia, PA, and approved May 7, 2007 (received for review February 16, 2007)
Abstract
We present improvements on our previous approaches for calculating vibrational spectroscopy observables for the OH stretch region of dilute HOD in liquid D_{2}O. These revised approaches are implemented to calculate IR and isotropic Raman spectra, using the SPC/E simulation model, and the results are in good agreement with experiment. We also calculate observables associated with threepulse IR echoes: the peak shift and 2DIR spectrum. The agreement with experiment for the former is improved over our previous calculations, but discrepancies between theory and experiment still exist. Using our proposed definition for hydrogen bonding in liquid water, we decompose the distribution of frequencies in the OH stretch region in terms of subensembles of HOD molecules with different local hydrogenbonding environments. Such a decomposition allows us to make the connection with experiments and calculations on water clusters and more generally to understand the extent of the relationship between transition frequency and local structure in the liquid.
Water is ubiquitous in science and nature (1), so it is natural that a tremendous amount of effort has been expended trying to describe and understand the structure and dynamics of its liquid state. Vibrational spectroscopy, both IR and Raman, provides an excellent probe of the local structure in water, because a local mode's vibrational frequency is exquisitely sensitive to the local mode's molecular environment. Actually, the cleanest information about local structure in water comes from the vibrational spectroscopy not of neat water, but rather of dilute HOD in either H_{2}O or D_{2}O, because in these cases, respectively, the OD or OH localmode stretch is almost completely decoupled from the other stretches in the liquid, thus functioning well as a local chromophore. IR and Raman spectra on these systems have been measured by many (2–9).
Valuable information about local dynamics in liquid water can also be obtained from vibrational spectroscopy experiments, in this case of the subpicosecond timedomain variety. On this time scale a local mode's vibrational frequency is continually changing because of molecular dynamics. The resulting dynamic frequency fluctuations, also known as spectral diffusion, can be measured by transient vibrational holeburning and threepulse echoes (5, 6, 10–22). In particular, these experiments provide information about both the shorttime (stretching) and longtime (making and breaking) aspects of intermolecular hydrogen bonds (23–35).
We and others have developed methods for the theoretical calculation of steadystate and ultrafast vibrational spectroscopy observables (7, 23–27, 36–40). In our approach, the single vibrational mode of interest, for example, the OH stretch of HOD (when it is immersed in D_{2}O), is treated quantum mechanically, whereas all other degrees of freedom (the bath) are treated classically. Thus we are making the adiabatic approximation in that for a given bath configuration one calculates the instantaneous eigenvalues of the OH stretch, leading to the relevant transition frequencies. The bath degrees of freedom evolve according to Newton's equations during a molecular dynamics simulation, generating a frequency trajectory from which one can calculate spectroscopic observables. To the extent that the transition dipole depends on the bath coordinates [nonCondon effects, which we think are important for water (33)] one also needs to obtain a transition dipole trajectory. The key, then, is to develop methods for calculating the transition frequency and dipole accurately and quickly for each bath configuration. For water we have suggested that both the frequency and dipole are approximately linearly related to the electric field from the surrounding waters, and the coefficients have been determined from quantum chemical calculations on clusters (25, 33, 35, 36). This approach has provided reasonable (although by no means perfect) agreement with experiment.
One of the purposes of this article is to report improvements on our previous methods. The improvements are 4fold: (i) our electronic structure benchmarks involve larger and more clusters; (ii) the electrostatic (point charge) effects of the water molecules surrounding each cluster are now included in the electronic structure calculations, ameliorating extrapolation issues; (iii) for each cluster anharmonic transition frequencies are calculated more accurately (using a discrete variable representation approach instead of fitting to a Morse potential); and (iv) a more accurate (quadratic instead of linear) fit of the benchmark frequencies to the electric field is proposed. With these improvements we recalculate experimental observables (Raman and IR spectra, threepulse echo peak shift, and 2DIR spectra) for HOD in D_{2}O by using the SPC/E water simulation model (41), finding in all cases modest improvements [over our previous calculations (35)] in the agreement with experiment.
The second purpose of this article is to understand better how the OH transition frequency is related to the Hbond configuration of the HOD molecule in the liquid. For years one has generally understood that molecules with weak or broken Hbonds (to the H of HOD) absorb on the blue side of the spectrum, whereas molecules with strong Hbonds absorb on the red. Renewed interest in a more detailed understanding of this general trend is motivated by the recent xray Raman, absorption, and emission experiments (42–49), some interpretations of which (46, 49) suggest a picture of Hbonding in liquid water quite different from the traditional one of more or less tetrahedral coordination with ≈3.5 Hbonds per molecule. The proposed picture is one of rings and chains of connected water molecules, each of which makes one strong donating and one strong accepting Hbond (46, 49). A crucial aspect of evaluating this scenario comes from assessing whether this picture is compatible with vibrational spectroscopy (50, 51). A second motivating factor for a better understanding of the relationship between transition frequency and local structure in the liquid comes from beautiful experiments and calculations on water clusters (52–59) showing clearly that vibrational frequencies of the OH stretch depend on much more than the geometry of the donor–acceptor pair; rather, they also depend on the number and nature of other Hbonds to the donor and acceptor molecules and indeed to the Hbonding configuration extending at least several solvation shells from the donor–acceptor pair. These complicated dependences arise from the cooperative nature of Hbonding interactions.
To shed some light on this issue, we decompose the distribution of OH stretch frequencies for the full ensemble of HOD molecules into subdistributions for HOD molecules with different Hbonding environments. Such a decomposition does not imply support for a mixture model of water; it simply is a way to classify instantaneous configurations. For such a decomposition one needs an Hbond definition, for which there are many proposals in the literature (46, 60–63). We have recently proposed an Hbonding definition (64) based on the electronic structure of water dimers (65), which involves the OH distance of the donor–acceptor pair, and the orientation of the acceptor molecule (in electronic structure parlance often called the electron donor). Thus we decompose the distribution of frequencies into subdistributions with different Hbonding configurations based on this definition, which has some relevance with regard to the compatibility of the IR and Raman spectra with recent xray absorption experiments and to the cluster results and their extrapolation to the liquid (57, 59).
OH Stretch Transition Frequency, Dipole, and Polarizability
In our approach to the vibrational spectroscopy of dilute chromophores in liquids, the first step is to generate clusters, in this case HOD molecules surrounded by D_{2}O molecules. To this end, we run a molecular dynamics simulation of 128 SPC/E D_{2}O molecules in the NVE ensemble, equilibrated to 300 K, as described (25, 33, 35, 36). Note that the experimental molecular number density at 1 atm pressure (1 atm = 101.3 kPa) is essentially the same for both water and heavy water; herein we use the value for heavy water of 3.32 × 10^{28}/m^{3}. For OH stretch frequency calculations we will imagine that one of the D atoms is an H. Thus we randomly tag a D atom on one of the molecules, and to generate clusters we retain all other water molecules whose O atoms are within 4 Å of the H (tagged D). Full electronic structure calculations will be performed on the molecules in the cluster. In addition, we record the positions of the nuclei of all other water molecules within half the box length; these nuclei will be point charges (keeping their SPC/E values) in the electronic structure calculations. Note, then, that this procedure for selecting clusters differs slightly from our previous one (one difference is that the average cluster size, now 9.25, is larger), and that in the past we did not include the point charges in the electronic structure calculations.
We then perform electronic structure calculations (B3LYP/6–311++G**) on 999 such randomly selected clusters, in the presence of the surrounding point charges. For these calculations all nuclei are fixed at their positions in the simulation, except that the OH stretch coordinate is varied (keeping the center of mass of the HOD molecule fixed). This maps out the 1D Born–Oppenheimer potential for the OH stretch, whose eigenvalues are solved by using the discrete variable representation and a pseudodiatomic HOD reduced mass of 0.954426 atomic mass units. Note that implementing this procedure for the isolated HOD molecule (OD distance and bond angle are kept at their SPC/E values), yields an OH stretch fundamental frequency of ω_{10} = 3,717 cm^{−1}, in good agreement with the experimental value of 3,707 cm^{−1}. In what follows, all calculated frequencies are then scaled by 0.9973 (obtained from the ratio of these two numbers) to correct for this small error.
As in the past, we now correlate these frequencies with the component of the electric field parallel to the OH bond, evaluated at the H atom, from the (SPC/E) point charges of all molecules in the configuration (except the HOD molecule) within half the box length. In the past we used a linear correlation. Here, however, we find that a quadratic dependence produces a slightly better fit. The results of this fit, compared with the ab initio results, are shown in Fig. 1, and the formula for the fit is given in Table 1. Although the fit is reasonable there still exists considerable scatter (rmsd is 70.1 cm^{−1}). Note that this rmsd is slightly larger than our earlier value (35), but now, because all molecules within half the box length are included in calculating the electric field, we have no residual extrapolation error as we did previously. We then rerun a simulation, randomly choosing D atoms as surrogate H atoms, calculate the electric field from all molecules within half the box length, and use the formula to relate field to frequency. This process generates a distribution of frequencies from the liquid, which is compared with the actual distribution of ab initio frequencies from the embedded cluster calculations (Fig. 2). The comparison is really very good, showing that even though there is quite of bit of inaccuracy from using the field to predict the frequency for individual configurations, on the average the formula produces a reasonably accurate distribution of frequencies. We also perform a similar quadratic fit for ω_{21}, which is needed for the nonlinear spectroscopy calculations (Table 1).
As in the past we write the projection of the transition dipole for vibrational states i and j in the direction of the external (light) electric field as: where μ′ is the dipole derivative, x is the OH stretch coordinate, û is the OH bond vector, and ϵ̂ is the direction of the external electric field. For the clusters μ′ is obtained from the electronic structure calculations, and the matrix elements x _{ij} are obtained from the discrete variable representation analysis. As before, μ′ and x _{ij} are then fit to the electric field and frequency ω_{ij}, respectively, with results shown in Table 1. These results are slightly different from those reported earlier, because of our somewhat different procedures for obtaining the frequencies. For the isotropic Raman spectrum we also need the polarizability; as before the polarizability derivative is obtained from electronic structure calculations and then correlated with the electric field as shown in Table 1.
IR, Raman, and 2DIR Spectra and Echo Peak Shift
Spectroscopic observables are calculated from frequency and transition dipole (polarizability) trajectories as in previous work (35). As such, nonCondon and motional narrowing effects are included. As before, the OH stretch fundamental lifetime is taken to be 700 fs (5). Experimental IR and Raman line shapes (5, 7), slightly corrected from the experimental crosssections and scattering intensities as described (35), are shown in Fig. 3, together with our theoretical calculations. Note that although the experimental Raman line shape is from the unpolarized signal (7) (whereas the theoretical result is for the isotropic line shape), in this particular instance the difference between these two is quite small (8, 9). The agreement between theory and experiment for the peak positions in both Raman and IR is excellent, hence the experimental red shift in going from Raman to IR is also properly described. The experimental difference in line shape between IR and Raman is also reproduced by the theory, although the theoretical shoulder in the Raman line shape is slightly too accentuated. The biggest discrepancy between theory and experiment is in the line widths, which are on the order of 10% too large in the theory. The comparison between theory and experiment shown herein represents a modest improvement over our previous results (35). A full discussion of the features of, and differences between, the IR and Raman line shapes has been given in ref. 35. Here, one might simply note that the substantial differences between the theoeretical line shapes and the frequency distribution (see Fig. 2) arise from nonCondon and motional narrowing effects.
Theoretical [calculated as described (35)] and experimental results for the threepulse echo peak shift (5) are shown in Fig. 4. One sees quantitative agreement with experiment at short (<50 fs) times and qualitative agreement thereafter. In particular, the theoretical recurrence occurs at 100 fs, whereas the experimental recurrence is at 150 fs (5), and the longtime decay of the theoretical peak shift is slightly too fast in comparison with experiment. Discrepancies between theory and experiment are caused by inadequacies of the SPC/E model and/or inaccuracies of our spectroscopic calculations. We are more inclined to believe the former, but certainly cannot rule out the latter. In any case, these theoretical peak shift results are again a modest improvement (in comparison with experiment) over our previous calculations for the same water model (35). Full 2DIR calculations for this model are shown in Fig. 5, which are qualitatively the same as those shown previously (35) and qualitatively the same as from experiment (20, 21). These 2DIR spectra can be analyzed in a number of different ways (66), including by considering the nodal slope between the 01 and 12 resonances (21, 35).
OH Stretch Frequency and HydrogenBonding Configuration
According to our theoretical calculations, the overall distribution of OH stretch frequencies in the liquid is given by the solid line in Fig. 2. As mentioned in the Introduction, we can decompose this distribution into subdistributions based on the Hbonding configurations of the HOD molecule. The classification of the configurations can be quite complex, involving several solvation shells of the HOD molecule. For this study, however, we will focus on the Hbonding state only of the HOD molecule itself. To this end, we can describe the Hbonding configuration of the HOD molecule by the number of Hbonds from neighboring molecules to the O atom, n _{O}, the number of Hbonds to the H atom, n _{H}, and the number of Hbonds to the D atom, n _{D}. Note that for configurations generated from an SPC/E simulation and reasonable Hbond definitions n _{O} = 0, 1, 2, 3 and n _{H} or n _{D} = 0, 1, 2. Thus every HOD molecule can be assigned to one of 36 classes.
This large number of classes makes it difficult to visualize and discuss results. However, because the probability of obtaining either n _{O} = 0 or 3 is small, we can profitably define a new number N _{O}, which is 1 if n _{O} is 0 or 1, and 2 if n _{O} is 2 or 3. Thus molecules with N _{O} = 1 correspond for the most part to those with one Hbond to the O, whereas molecules with N _{O} = 2 correspond for the most part to those with two Hbonds to the O. Likewise, the probability of obtaining n _{H} = 2 is small, and so we can define a new number N _{H}, which is 0 if n _{H} = 0, and 1 if n _{H} = 1 or 2, and similarly for D. This then gives a more manageable number of eight Hbond classes, each of which is described by the triplet of numbers N _{O}, N _{H}, and N _{D}.
Alternatively, we can describe each class by the total number of Hbonds N _{O} + N _{H} + N _{D}, the number of Hbond donors N _{H} + N _{D}, and if there is a single donor (in which case the latter number is 1), whether it is the H or the D. Thus, for example, the triplet N _{O} = 2, N _{H} = 1, and N _{D} = 0 can be labeled 3_{SH} (three Hbonds with the single donor being the H). (Note, however, that because of the way we have combined the original 36 classes, a member of the class 3_{SH} will occasionally have four Hbonds.) The translation between these two different ways of labeling environments is given in Table 2, where D (as the first letter) means double donor, S means single donor, and N means nondonor.
To perform this classification for an HOD molecule in the liquid, we need to be able to specify how many Hbonds each of its atoms is involved with. For this we turn to our proposed Hbond definition (64) based on the electronic occupancy of OH antibonding orbitals (67). For a water monomer the occupancies of its two OH antibonding orbitals are essentially zero. For an H atom (donor) involved in an Hbond, however, the occupancy of its OH antibonding orbital can be appreciable, because of electron donation from the lone pairs of the acceptor molecule (67). For a given pair of H and O atoms on different molecules we have shown, using electronic structure calculations, that this occupancy depends primarily on two variables related to the relative geometry of these two molecules (64). One variable is r, the HO intermolecular distance. The second is the angle that the intermolecular OH ray makes with the outof plane unit vector (in the direction of the p _{z} orbital) on the acceptor molecule (64). ψ is defined to be the smaller of this angle and its complement. The map is given by (64): where N is the occupancy and ψ is in degrees. We have calculated (64) the distribution of occupancies for all intermolecular OH pairs in the liquid, and from its bimodal structure we have suggested a critical occupancy of 0.0085. Therefore, for every (intermolecular) pair of H and O atoms in the H_{2}O liquid, one defines the coordinates r and ψ and calculates the occupancy N according to the above; if N > 0.0085, these atoms are Hbonded (and if N < 0.0085, they are not). According to this definition, for liquid H_{2}O at 300 K and 1 atm pressure 〈n _{H}〉 = 0.84, and 〈n _{O}〉 = 2 〈n _{H}〉 = 1.68, so that the average number of Hbonds associated with a given molecule is 〈n _{O}〉 + 2 〈n _{H}〉 = 3.36 (64). Furthermore, the fractions of molecules with one, two, three, four, or five Hbonds are 0.017, 0.126, 0.372, 0.449, and 0.036, respectively.
Finally, then, for the purpose of the frequency calculations for HOD in D_{2}O, we randomly tag D atoms in the D_{2}O simulation to be surrogate H atoms. For a given configuration of the HOD molecule in the liquid, N _{O}, N _{H}, and N _{D} are determined as described above. At the same time the OH stretch frequency is calculated as described above. This process is repeated many times, leading to a distribution of OH stretch frequencies for each Hbonded class. The eight subdistributions are shown in Fig. 6, together with the overall frequency distribution. The subdistributions are not normalized individually; thus subdistributions sum to give the overall distribution. Each subdistribution can be characterized by its average frequency, and the fraction of molecules it contains, as shown in Table 2.
From Fig. 6 or Table 2 one can see what happens to the OH stretch frequency when a hydrogen bond to the H is added (N _{H} goes from 0 to 1, keeping N _{O} and N _{D} unchanged): the average OH stretch frequency is redshifted by anywhere from 192 to 229 cm^{−1}. Similarly, when an Hbond to the O is added (N _{O} goes from 1 to 2), the average frequency is redshifted by anywhere from 41 to 77 cm^{−1}, and when an Hbond to the D is added (N _{D} goes from 0 to 1), the average frequency is blueshifted by anywhere from 11 to 19 cm^{−1}. Consequently it follows that the 2_{SD} distribution is the least redshifted (from the gasphase frequency of 3,707 cm^{−1}), and 3_{SH} is the most redshifted. The frequency ordering of all of the distributions can be read off of Table 2. Interestingly, this ordering shares some features with that from the recent paper by Lenz and Ojamäe (57) obtained with an extrapolation from calculations on H_{2}O clusters. In particular, the frequency ordering 3_{SH} < 4_{D} < 3_{D} is found in both cases. The same conclusion is also found in the interesting theoretical and experimental paper on water clusters by Ohno et al. (59). In light of the results in both of these cluster papers it seems likely that the large breadth of each of our four Hbonded (to the H) subensembles can be attributed in part to different Hbonding configurations of the acceptor molecule.
It is important to reiterate that for each HOD molecule the frequency calculation is obtained from the electric field fit, which has significant scatter compared with the benchmark frequencies, as discussed above. Thus it would be better if we could use the benchmark frequencies themselves (rather than those from the electric field map) to investigate these distributions. Evidently, however (see Fig. 2), 999 clusters are not enough to obtain even a good overall distribution of frequencies, much less good subdistributions. To obtain reasonable subdistributions from benchmarks would take tens of thousands of clusters. Also note that, in general, line shapes are not simply frequency distributions, because of the complications of nonCondon effects and motional narrowing. Despite these realities, we believe that one can draw two further conclusions from the above. First, from summing the relative probabilities of the different classes, we find that the fraction of double donors is 0.702, the fraction of single donors is 0.272, and the fraction of nondonors is 0.027. If the fraction of single donors were substantially higher, as claimed by Wernet et al. (46), and of course making the big assumption that the shapes and positions of these subdistributions are invariant to the particular microscopic model of water, the amplitudes of the four single donor distributions would be increased at the expense of the two double donor distributions, leading to a dramatic change in the shape of the overall distribution. In particular, the 2_{SD} and 3_{SD} distributions peak at the relatively high values of 3,665 and 3,624 cm^{−1}, respectively, because these correspond to “dangling” OH groups. A large fraction of single donors would then imply strong spectral intensity ≈3,650 cm^{−1}, and maybe even a doubly peaked spectrum (51), which is not seen experimentally. Second, the eight subdistributions are broad and overlapping, implying that the structural origins of the line shapes are complicated. A corollary, of course, is that detailed structural inferences from line shape features must be drawn with care.
Summary
Advances in our theoretical approach for describing liquidstate spectroscopy have led to improved results (in comparison with experiment) for 1D and 2D line shapes for dilute HOD in liquid D_{2}O. The theoretical distribution of frequencies has been decomposed into subensembles for different Hbonding environments of the HOD molecules, each of which peak at a different frequency. The still quite substantial breadth of four of these subdistributions is presumably caused in part by different Hbonding configurations of the acceptor and other D_{2}O molecules. Our theoretical results indicate that the fraction of HOD molecules that are single donors is ≈27%, and we suggest that if, in fact, this fraction was significantly higher, the resulting vibrational line shapes would not be compatible with experiment.
Acknowledgments
We thank Prof. Will Polik and Hope College (Holland, MI) for the use of their MU3C computer cluster on which some of the calculations reported herein were performed and Prof. Lars G. M. Pettersson and Mathias Ljungberg for helpful discussion. This work was supported by National Science Foundation Grant CHE0446666, an American Chemical Society/Petroleum Research Fund AC grant, and a Hertz fellowship (to J.R.S.).
Footnotes
 ^{†}To whom correspondence should be addressed. Email: skinner{at}chem.wisc.edu

Author contributions: J.L.S. designed research; and B.A., R.K., and J.R.S. performed research.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.
 © 2007 by The National Academy of Sciences of the USA
References

↵
 Ball P

↵
 Falk M ,
 Ford TA
 ↵
 ↵
 ↵
 ↵

↵
 Smith JD ,
 Cappa CD ,
 Wilson KR ,
 Cohen RC ,
 Geissler PL ,
 Saykally RJ
 ↵
 ↵
 ↵

↵
 Woutersen S ,
 Emmerichs U ,
 Bakker HJ
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Fecko CJ ,
 Eaves JD ,
 Loparo JJ ,
 Tokmakoff A ,
 Geissler PL
 ↵

↵
 Eaves JD ,
 Loparo JJ ,
 Fecko CJ ,
 Roberts ST ,
 Tokmakoff A ,
 Geissler PL
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Schmidt JR ,
 Roberts ST ,
 Loparo JJ ,
 Tokmakoff A ,
 Fayer MD ,
 Skinner JL
 ↵
 ↵
 ↵
 ↵

↵
 Harder E ,
 Eaves JD ,
 Tokmakoff A ,
 Berne BJ
 ↵

↵
 Wilson KR ,
 Rude BS ,
 Catalano T ,
 Schaller RD ,
 Tobin JG ,
 Co DT ,
 Saykally RJ
 ↵
 ↵
 ↵

↵
 Wernet P ,
 Nordlund D ,
 Bergmann U ,
 Cavalleri M ,
 Odelius M ,
 Ogasawara H ,
 Näslund LÅ ,
 Hirsch TK ,
 Ojamäe L ,
 Glatzel P ,
 et al.

↵
 Smith JD ,
 Cappa CD ,
 Wilson KR ,
 Messer BM ,
 Cohen RC ,
 Saykally RJ
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Lenz A ,
 Ojamäe L
 ↵
 ↵
 ↵
 ↵
 ↵

↵
 Kuo IFW ,
 Mundy CJ
 ↵
 ↵
 ↵

↵
 Weinhold F ,
 Landis C