16299
2009, 113, 16299–16302 Published on Web 08/26/2009
Vibrational Spectroscopic Response of Intermolecular Orientational Correlation at the Water Surface Tatsuya Ishiyama and Akihiro Morita* Department of Chemistry, Graduate School of Science, Tohoku UniVersity, Sendai 980-8578, Japan ReceiVed: June 30, 2009; ReVised Manuscript ReceiVed: August 14, 2009
The second-order nonlinear susceptibility, which is the source of the sum frequency generation (SFG) spectrum, for the vapor-water interface is calculated by molecular dynamics simulation based on the time correlation formula. The calculated susceptibility well reproduces that of the recent phase-sensitive SFG experiment (Ji et al. Phys. ReV. Lett. 2008, 100, 096102), including the positive imaginary part in the 3000-3200 cm-1 region, which is a controversial problem in the assignment. It is found that this positive part of the imaginary susceptibility is induced by the intermolecular orientational correlation at the interface through the anisotropic local field. Water interfaces play crucial roles in many fields of science, including atmospheric physical and chemical processes, bubble engineering, protein stabilities and functions, etc. Because of its important but intricate nature, a large number of studies about water interfacial structure have been carried out both experimentally1-6 and theoretically.7-14 Vibrational sum frequency generation (SFG) spectroscopy is a powerful probe of the hydrogen-bond structure and molecular orientation at the water interface.15,16 Regarding the orientational structure of surface water, one of the most noteworthy topics in recent studies is the phase-sensitive measurements of SFG,17-21 which can provide phase information of the nonlinear susceptibility χ, whereas usually the observed intensity of SFG spectra is governed by |χ|2. The phase of χ, in particular the sign of the imaginary part of χ, Im[χ], is considered to reflect the polar orientation of surface water.22 The phase-sensitive SFG measurement of the vapor-water interface was first reported by Ji et al. in 2008,18 revealing three definite components in the imaginary part of χ (as discussed later in Figure 2a). The first is the positive region at about 3700 cm-1, which is assigned to the dangling O-H moieties that direct their dipoles toward the vapor phase. The second is the negative region in the 3200-3600 cm-1 range, which is attributed to the hydrogen-bonded O-H moieties that direct their dipoles toward bulk liquid. The third is the positive region in the 3000-3200 cm-1 range. The third positive region has been surprising, since none of the theoretical calculations conducted so far have been able to predict this component.7-14 Elucidating this positive region is of key importance in the recent controversy in the assignment of the water SFG spectrum,23,24 and Tian and Shen requested theorists to try to simulate this positive region.20 In the present Letter, this puzzling issue is clarified by appropriately incorporating the local field effect in the molecular dynamics (MD) simulation, that the anisotropic local field renders Im[χ] positive in the 3000-3200 cm-1 region. The mechanism of the anisotropic * To whom correspondence should be addressed. E-mail: amorita@ mail.tains.tohoku.ac.jp.
10.1021/jp9060957 CCC: $40.75
local field is explained from the viewpoint of intermolecular orientational correlation. In the following, the improved treatment of the local field in the SFG calculation is formulated and validated. The second-order nonlinear susceptibility χ ()χR + χNR) consists of the vibrational resonant term χR and the nonresonant term χNR. The former, χR, is given using a time correlation function of the polarizability tensor Apq and the dipole vector Mr of the interface system as follows:12 R χpqr )
iωIR kBT
∫0∞ dt exp(iωIRt)〈Apq(t)Mr(0)〉
(1)
where kB and T are the Boltzmann constant and temperature and 〈〉 denotes the statistical average. A and M are calculated by polarizable MD simulation, where electric field and polarization are determined self-consistently at each time step. The nonresonant term χNR, on the other hand, is treated as a constant over the range of infrared frequency ωIR in question. Throughout this paper, the suffixes p, q, r, and s refer to the Cartesian components (x-z) and i and j to molecules. Suppose the ith water molecule has a polarizability Ri and dipole moment (permanent + induced) µi. Then, the selfconsistent condition between the dipole moment µi and the 0 + electric field at the ith molecule Ei is given as µp,i ) µp,i 0 - ∑j(*i) ∑qx,y,z Tp,i;q,jµq,j, where µ0 is ∑qx,y,z Rpq,iEq,i and Ep,i ) Ep,i the permanent dipole moment, E0i is the electric field generated by the surrounding permanent charges, and Tij is the dipole-dipole interaction tensor.25 With the interaction site point charge model, the electric field E0i is explicitly written as sites
Ei0 )
∑ ∑
j(*i)
f(rija)Qja
a
rija3
rija
(2)
where Qja denotes the partial charge on site a of the jth molecule located at rja, rija ) ri - rja, and rija ) |rija|. f(r) is the damping 2009 American Chemical Society
16300
J. Phys. Chem. C, Vol. 113, No. 37, 2009
Letters
Figure 1. (b) Dipole moment and (c) potential energy surface of the water dimer illustrated in panel a, where the direction of the x axis is defined along the two centers of mass of the water molecules. In panels b and c, the solid lines with filled circles stand for the DFT calculations, while the blue and red symbols stand for the polarizable molecular model with f ) 1 and f ) f Thole, respectively.
factor as a function of distance r. The damping factor f(rija) e 1 in eq 2 reduces the electric field generated by the point partial charge Qja in a certain short distance, and it will be focused on below. From the above conditions, the dipole moment M of the whole interface system can be expressed as25 x,y,z
Mp )
∑ ∑ Fr,piµr,i′ i
(3)
r
0 0 where µr,i ′ ) µr,i + ∑q Rrq,iEq,i . Fr,pi is the local field correction i factor formulated as Fr,p ) ∑j [G-1]ri,pj and [G]ri,pj ) δprδij + ∑q Tr,i;q,jRqp,j with δij being the Kronecker’s delta. The system polarizability Apq in eq 1 is also expressed by Apq ) ∑j ∑r ∑s Fr,pjFs,qjRrs,j. The damping factor f in eq 2 is necessary for the classical polarizable MD simulation to be stable. If one employs f ) 1 (no damping), the MD trajectories often suffer from unphysical divergence of polarization, called “polarization catastrophe”.26 A number of polarizable force fields, including our point dipole (PD) model,27 have employed the Thole-type damping function f ) f Thole(r)26 or its analogues. While it is generally believed that the damping factor f remedies the shortcomings of shortrange Coulomb interaction by smearing the charge distribution of the point charge/dipole model, its physical meaning or its quantitative accuracy is not clear enough. We therefore examine the damping factor in describing the short-range interaction in comparison with density functional theory (DFT). The potential energy and the dipole moment of the water dimer with varying water-water distance are predicted by the PD water model27 with two damping functions, f ) f Thole(r) and f ) 1, and those results are compared to B3LYP/d-aug-cc-pVDZ28,29 calculations using Gaussian 03.30
Shown in Figure 1a is the equilibrium configuration of the water dimer calculated by DFT, where the O-H distance of the hydrogen bond is rOH ) 1.946 Å. In Figure 1b, the dipole moments of the dimer along the hydrogen-bonding O-H direction (see Figure 1a and its caption) are shown as a function of rOH, by moving the donor hydrogen with the other nuclei fixed. We note that Figure 1b is relevant to the transition moment for the O-H stretching vibration under the hydrogenbonding environment. One can see that the case of f ) 1 well reproduces the dipole moment in the DFT calculation, while the case of f ) f Thole describes poor Mx. This demonstrates that the damping factor is not necessary to well reproduce the electrostatic properties such as the dipole moment, as the original point charge/dipole model (without damping)27 has been optimized to describe the electrostatic field around the molecule. On the other hand, Figure 1c shows the intermolecular potential energy surface plotted with varying rOH, where the molecular geometry of each water is fixed. This panel clearly shows that the case of f ) f Thole agrees well with the DFT calculations. The result suggests that the appropriate force field for intermolecular interaction should be the case of f ) f Thole. In summary, while the force field for the polarizable MD simulation should use f ) f Thole, the perturbation in electrostatic properties by strong hydrogen bond can be described more appropriately with f ) 1. This is an important point that the previous theoretical calculations have failed to notice. The mechanism of this difference will be discussed in detail elsewhere.31 Next, we calculated M and A during the MD simulation of the water surface27,32 with two damping factors, f ) f Thole and f ) 1, for comparison, while in both cases the MD trajectories were generated using f ) f Thole. Then, χR is derived with eq 1 in both cases, while the nonresonant term χNR of the water surface is assumed a real constant value over the range of infrared frequency, so as to be consistent with the experimental observation.33 Shown in Figure 2 is the frequency-dependent nonlinear susceptibility χqqz for the water surface, where q stands for x or y tangential to the surface and z is the surface normal direction from the liquid to vapor phase. Panel a is the experimental result by Ji et al.,18 whereas panels b and c show the calculated results with f ) f Thole and f ) 1, respectively. One immediately notices that panel c (f ) 1) rather than panel b (f ) f Thole) actually reproduces the positive region of Im[χqqz] in the 3000-3200 cm-1 range. This result also supports the validity of no damping treatment (f ) 1) for the calculation of the system electrostatic properties. In what follows, we elucidate that the positive Im[χqqz] in the 3000-3200 cm-1 range originates from the anisotropic local field effect. For this purpose, the system dipole moment of eq 3 in the z direction, Mz, is decomposed into the isotropic part Mziso and the anisotropic part Mzaniso: Mz ) Mziso + Mzaniso, where
Mziso )
∑ Fz,ziµz,i′ , i
Mzaniso )
∑ (Fx,ziµx,i′ + Fy,ziµy,i′ ) i
(4) aniso Accordingly, χqqz is decomposed into χqqz ) χiso qqz + χqqz , where ∞ iso aniso ∞ ∼ ∫ dt exp(iω t)〈A (t)M (0)〉 and χ ∼ ∫ χiso dt exp(iωIRt)qqz 0 IR qq z qqz 0 〈Aqq(t)Mzaniso(0)〉. The imaginary parts of two components, iso aniso ] and Im[χqqz Im[χqqz ], are shown in Figure 3. One can see that iso ] (shown in a blue line) the isotropic component Im[χqqz essentially consists of the positive region above 3600 cm-1 and the negative region below 3600 cm-1, which reflects the averaged direction of the dangling and hydrogen-bonded OH moieties, respectively, near the interface. An interesting result
Letters
J. Phys. Chem. C, Vol. 113, No. 37, 2009 16301
Figure 4. Average dipole moment in the vicinity of the water interface (see the text). zˆ ) 0 is the Gibbs dividing surface, zˆ < 0 is the liquid side, and zˆ > 0 is the vapor side: (a) isotropic part; (b) anisotropic part.
Figure 2. Real (blue) and imaginary (red) parts of χqqz in the ssp polarized SFG spectrum for the water surface: (a) experimental results by Ji et al.;18 (Copyright 2008, American Physical Society); (b) calculated results with f ) f Thole; (c) calculated results with f ) 1.
Figure 3. Decomposed result of Im[χqqz] calculated with f ) 1 into the isotropic component Im[χiso qqz] (blue line) and anisotropic component aniso ] (red line). Im[χqqz aniso is that the red line, the anisotropic component Im[χqqz ], is -1 positive from 3000 to 3400 cm , which makes the total Im[χqqz] positive in the 3000-3200 cm-1 region. We further investigated the spatial origin of each component, Mziso or Mzaniso in eq 4, by mapping the suffix i into the depth zˆ where the molecule i is located. Figure 4 shows the average dipole moments of isotropic ′ 〉 and 〈Fx,ziµx,i ′ + Fy,ziµy,i ′ 〉 in and anisotropic components, 〈Fz,ziµz,i eq 4, as a function of the interfacial depth zˆ, demonstrating that the isotropic component (panel a) is negative near the Gibbs dividing surface (zˆ ) 0),34 while the anisotropic component (panel b) is positive there. This result indicates that the anisotropic correlation is present only near the top monolayer below the Gibbs dividing surface. The mechanism of the anisotropic component is schematically ′ >0 illustrated in Figure 5. The positive Mzaniso implies Fx,ziµx,i ′ for a molecule i at the interface, indicating that Fx,zi > 0 and µx,i
Figure 5. Schematics of the anisotropic local field and intermolecular correlation: (a) condition of Fx,zi > 0 and µx,i ′ > 0 for the surface water; (b) orientational correlation of adjacent water molecules with Rxz,j > 0, generating the anisotropic local field (see the text).
> 0 or vice versa. (Note that the case of Fx,zi < 0 and µx,i ′ < 0 is equivalent if one reverses the x axis.) The question addressed ′ tend to have the same sign. The here is why Fx,zi and µx,i anisotropic local field factor of Fx,zi > 0 means that when an external field along z is imposed (Ez0 > 0), then the local field along x is induced (Ex > 0). This anisotropy originates from the anisotropic orientation of the surrounding water molecules. ′ > 0), it When the molecular dipole of i is laterally oriented (µx,i tends to induce a lateral alignment of the adjacent dipoles through the hydrogen-bond network. This local alignment of the water orientation induces positive Rxz,j, as illustrated in panel b, which creates the anisotropic local field of Fx,zi > 0, and hence positive Mzaniso. We confirmed this mechanism by investigating surface structure by the MD simulation, and found that the intermolecular orientational correlation along the x direction becomes significant in a strongly hydrogen-bonding pair of adjacent water molecules. Such strong hydrogen-bonding pairs make the above mechanism of anisotropic local field particularly significant in the strongly red-shifted region for 3000-3400 cm-1. In the present study, we calculated the frequency-dependent second-order susceptibility χ of the water surface by MD simulation, and elucidated the phase behavior of χqqz reported by the recent experiment. In particular, the positive Im[χqqz] in the 3000-3200 cm-1 range was reproduced, which is a recent controversial problem in the assignment. Our MD analysis revealed that this positive imaginary susceptibility is not necessarily attributed to the orientation of individual molecular dipoles pointing to the vapor phase but originates from intermolecular orientational correlation of surface water, which results in the anisotropic local field as demonstrated above. This anisotropic local field effect is analogous with the previously
16302
J. Phys. Chem. C, Vol. 113, No. 37, 2009
reported intermolecular vibrational correlation effect,35 in a sense that the strong dipole-dipole coupling of aqueous interfaces makes the spectral interpretation in terms of collective motions. This strong coupling can be properly described by using the accurate damping factor f ) 1 for the calculation of M and A, and this damping factor is validated in comparison with DFT calculations. The present study not only solves the recent controversial problem on the phase of χ concerning the water surface but also indicates that the interpretation on such SFG spectra should be refined in terms of collective behavior of molecules. Acknowledgment. We are grateful to Dr. Tahei Tahara for valuable comments. This work was supported by the NextGeneration Supercomputer Project and the Grants-in-Aid (Nos. 20750003, 20038003, 20050003, and 21245002) by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. References and Notes (1) Du, Q.; Superfine, R.; Freysz, E.; Shen, Y. R. Phys. ReV. Lett. 1993, 70, 2313. (2) Du, Q.; Freysz, E.; Shen, Y. R. Science 1994, 264, 826. (3) Shultz, M. J.; Schnitzer, C.; Simonelli, D.; Baldelli, S. Int. ReV. Phys. Chem. 2000, 19, 123. (4) Raymond, E. A.; Tarbuck, T. L.; Brown, M. G.; Richmond, G. L. J. Phys. Chem. B 2003, 107, 546. (5) Liu, D.; Ma, G.; Levering, L. M.; Allen, H. C. J. Phys. Chem. B 2004, 108, 2252. (6) Sovago, M.; Campen, R. K.; Wurpel, G. W. H.; Mu¨ller, M.; Bakker, H. J.; Bonn, M. Phys. ReV. Lett. 2008, 100, 173901. (7) Benjamin, I. Phys. ReV. Lett. 1994, 73, 2083. (8) Buch, V. J. Phys. Chem. B 2005, 109, 17771. (9) Brown, E. C.; Mucha, M.; Jungwirth, P.; Tobias, D. J. J. Phys. Chem. B 2005, 109, 7934. (10) Perry, A.; Neipert, C.; Space, B.; Moore, P. B. Chem. ReV. 2006, 106, 1234. (11) Walker, D. S.; Hore, D. K.; Richmond, G. L. J. Phys. Chem. B 2006, 110, 20451. (12) Morita, A.; Ishiyama, T. Phys. Chem. Chem. Phys. 2008, 10, 5801. (13) Auer, B. M.; Skinner, J. L. J. Phys. Chem. B 2009, 113, 4125.
Letters (14) Noah-Vanhoucke, J.; Smith, J. D.; Geissler, P. L. J. Phys. Chem. B 2009, 113, 4065. (15) Richmond, G. L. Chem. ReV. 2002, 102, 2693. (16) Shen, Y. R.; Ostroverkhov, V. Chem. ReV. 2006, 106, 1140. (17) Ostroverkhov, V.; Waychunas, G. A.; Shen, Y. R. Phys. ReV. Lett. 2005, 94, 046102. (18) Ji, N.; Ostroverkhov, V.; Tian, C. S.; Shen, Y. R. Phys. ReV. Lett. 2008, 100, 096102. (19) Tian, C. S.; Shen, Y. R. Chem. Phys. Lett. 2009, 470, 1. (20) Tian, C. S.; Shen, Y. R. J. Am. Chem. Soc. 2009, 131, 2790. (21) Nihonyanagi, S.; Yamaguchi, S.; Tahara, T. J. Chem. Phys. 2009, 130, 204704. (22) Morita, A.; Hynes, J. T. Chem. Phys. 2000, 258, 371. (23) Tian, C. S.; Shen, Y. R. Phys. ReV. Lett. 2008, 101, 139401. (24) Sovago, M.; Campen, R. K.; Wurpel, G. W. H.; Mu¨ller, M.; Bakker, H. J.; Bonn, M. Phys. ReV. Lett. 2008, 101, 139402. (25) Morita, A.; Hynes, J. T. J. Phys. Chem. B 2002, 106, 673. (26) Thole, B. T. Chem. Phys. 1981, 59, 341. (27) Ishiyama, T.; Morita, A. J. Phys. Chem. C 2007, 111, 721. (28) Dunning, T. H., Jr. J. Chem. Phys. 1989, 90, 1007. (29) Kendall, R. A.; Dunning, T. H., Jr. J. Chem. Phys. 1992, 96, 6796. (30) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision D.02; Gaussian, Inc.: Wallingford CT, 2004. (31) Ishiyama, T.; Morita, A. To be submitted. (32) Ishiyama, T.; Morita, A. J. Phys. Chem. C 2007, 111, 738. (33) Morita, A. J. Phys. Chem. B 2006, 110, 3158. (34) Adamson, A. W.; Gast, A. P. Physical Chemistry of Surfaces, 6th ed.; Wiley: New York, 1997. (35) Ishiyama, T.; Morita, A. Chem. Phys. Lett. 2006, 431, 78.
JP9060957
