Effects of hydrophobic solute on water normal modes
Abstract
The vibrational modes of water get significantly modified by external solutes; this becomes particularly important when the solute is hydrophobic. In this work, we examine the effects of a tiny hydrophobe, methane, on the normal modes of water, using small cluster-based harmonic normal mode analysis of aqueous methane system. We estimate the vibrational density of states and also the infrared spectral density. We compare the methane-water data with the bulk water response. We decompose these modes based on different vibrational characters. The stretch-bend decomposition reflects a pronounced coupling between the methane asymmetric stretch mode and the water symmetric stretch mode. We examine the methane-water data in terms of the symmetry of the central water molecule’s vibrations and find that asymmetric modes do not contribute. We also find that the vibrational modes having non-zero contributions from methane molecule are extremely localized in nature.
I Introduction
Hydrophobic hydration plays a fundamental role in several biological processes such as membrane formation, protein folding, ligand binding, assembly of proteins into functional complexes, lipid bi-layer formation, and many others Tanford1973 ; Srinivas2002 ; Srinivas2003 ; Mason2004 ; Chandler2005 ; Dyson2006 ; Bakulin2011 ; Hazra2014 ; Ben2015 ; Grdadolnik2017 . It is also evident that the abundant hydrogen bond (H-bond) network of bulk water gets significantly modified in the presence of different solutes, particularly important in the case of hydrophobic molecules. The influence of such molecules on the structure of liquid water is an essential topic of ongoing research. The quest for hydrophobic hydration has emerged with the seminal work of Frank and Evans Frank1945 who anticipated that “water forms frozen patches or microscopic icebergs around such solute molecules.” The low solubility of any such oily solutes poses significant experimental challenges to probe the system experimentally. Thus a numerous earlier explorations of hydrophobic hydration in the case of the classic hydrophobe, methane (), have been done mainly by employing theoretical,Kauzmann1959 molecular dynamics (MD) Rossky1982 ; Laza1992 ; Ash1996 ; Sharp1997 ; Hummer1998 ; Pas2009 , and ab initio (AIMD) simulation techniques Mon2012 ; Galamba2014 ; Grdadolnik2017 .
A few experimental findings of hydrophobic hydration portray conflicting conclusions Dejong1997 ; Buchanan2005 ; Laage2009 ; Perera2009 ; Tie2010 ; Str2014 ; Rankin2015 ; Grdadolnik2017 ; Ben2018 . The neutron diffraction study by De Jonget al. concluded that 19 tangentially oriented water molecules surround the methane molecule Dejong1997 , and a similar experiment by Buchananet al. concluded “no evidence for enhanced water structure.” Buchanan2005 A recent work by Grdadolnik et al. using infrared (IR) spectroscopy of isotopically dilute systems in addition to ab initio MD suggested “evidence that supports, unequivocally, the iceberg view of hydrophobicity.” This is inconsistent with other classical MD Rossky1982 ; Laza1992 ; Ash1996 ; Sharp1997 ; Hummer1998 ; Pas2009 and AIMD Mon2012 ; Galamba2014 ; Grdadolnik2017 data which find no icebergs around hydrophobic molecules. However, the MD hydration thermodynamic analyses Qvist2008 exactly reproduce the entropic anomalies that originally encouraged the Frank and Evans iceberg hypothesis Frank1945 . Furthermore, experimental fs-IR Rezus2007 ; Petersen2009 ; Str2014 ; Rankin2015 ; Ben2018 , NMR Laage2009 , and theoretical explorations Ros2012 ; Mon2012 find slower water dynamics around hydrophobic moieties, with no evidence of iceberg formation. Recent Raman multivariate curve resolution (Raman-MCR) hydration-shell vibrational spectroscopic measurements find that “methane’s hydration-shell is slightly more clathrate-like than pure water at ambient, and lower temperatures.” Ben2018 ; Ben2019 ; Bredt2020
Despite the previous conflicting conclusions, it is evident that the underlying many-body interactions could give rise to different hydration arrangements around the hydrophobe Wiggins1997 ; Lum1999 ; Huang2000 ; Chandler2005 ; Bakulin2011 ; Galamba2013 ; Galamba2014 ; Grdadolnik2017 . This modifies the strength of the hydrogen bonds in the close vicinity of the hydrophobe Scatena2001 ; Raschke2005 ; Buchanan2005 ; Bakulin2011 ; Davis2012 ; Mon2012 ; Grdadolnik2017 . In one of the earlier works, Stillinger has argued that for sufficiently weak solute-water attraction, a large smooth hydrophobic surface might be enclosed by a microscopically thin film of water vapor Stillinger1973 .
The high sensitivity of the IR spectroscopy to local structures makes it an appealing method for exploring of these systems, where the hydrophobic moieties alter the water H-bond significantly Fecko2003 ; Devendra2020 . IR spectroscopy of O–H stretching mode serves as the most reliable and sensitive probe for investigating the relative strengths of hydrogen bonds Davis2012 ; Fecko2003 ; Hecht1992 ; Hecht1993 ; Sharp1997 ; Auer2007 ; Laage2009 ; Stiopkin2011 ; Biswas2013 ; Biswas2016 . Despite its extensive use, the superposition of several solvation structures embedded in the experimental data makes it extremely difficult to establish the structure-spectrum correlation accurately. Besides, a further delocalization arises due to the strong inter molecular an-harmonic couplings. In contrast, the molecular level resolution of computer simulation aided theoretical spectroscopic techniques enables us to probe these systems at the microscopic level Auer2007 ; Biswas2013 ; Biswas2016 ; Corcelli2004 ; Corcelli2005 ; Roberts2009 ; Biswas2017 ; Samanta2018 . In several earlier works, mixed quantum-classical (MQC) spectroscopy models have been extensively used to study the isotope dilute aqueous systems using O–H stretching vibrations as probe Auer2007 ; Biswas2013 ; Biswas2016 ; Corcelli2004 ; Corcelli2005 ; Samanta2018 . However, this reduced dimensional description does not allow us to investigate the extent of mode mixing in the condensed phase. The instantaneous normal mode analysis is found to be extremely useful in this regard Biswas2017 ; Cho1994 . In this work, we explore the methane-water system’s vibrational responses using instantaneous normal mode analyses of small methane-water clusters drawn from classical MD trajectories. A particular emphasis is given to understand the effects of the hydrophobic molecule on the mode-mixing and the extent of delocalization.
The rest of this paper is organized as follows: we discuss the simulation details, cluster selection, and normal mode calculation details in Section II, Section III includes the results and discussions, and conclusions are given in Section IV.
II Simulation Details
We use the Groningen machine for chemical simulations (GROMACS) version 2019.1 Abraham2015 to perform classical molecular dynamics simulations of two systems: (1) methane-water system consisting of one methane molecule and 255 water molecules, and (2) bulk water system with 256 water molecules. We use optimized potentials for liquid simulations all-atom (OPLS-AA) Jorgensen1996 force field to describe the methane molecule, and extended simple point charge (SPC/E) Berendsen1987 force field to describe the water molecules. We employ periodic boundary conditions in all three directions. First, we perform energy minimization of the systems using the steepest descent algorithm, followed by 1 ns equilibration in the NPT ensemble. Thereafter, the 5 ns production run is carried out in the NVT ensemble. We use Berendsen thermostat Berendsen1984 with a relaxation time of 0.1 ps, and Parrinello–Rahman barostat Parrinello1981 ; Nose1983 with a relaxation time of 1.0 ps for maintaining a constant temperature at 300 K, and pressure at 1 bar, respectively. For neighbour searching, and non-bonded interactions, we use a 9 Å cut-off radius, and all the bonds are kept fixed using linear constraint solver (LINCS) Hess1997 algorithm. We use Particle Mesh Ewald (PME) Darden1993 with an FFT grid spacing of 1.6 Å for the long-range electrostatic interaction’s calculation.

We perform the harmonic normal mode analysis for both systems by employing the following strategy. We utilize electronic structure calculations on instantaneous small frozen clusters extracted from the classical MD trajectory. From the classical trajectory, we extract small clusters by identifying a central O atom that belongs to the water molecule closest to the methane molecule and including all water molecules having their oxygen atoms within a 5.0 Å radius (Figure 1). We obtain the methane-water data by using such 2000 methane-water clusters. We choose 1500 small clusters extracted similarly from the bulk water trajectory to extract the bulk water response. The methane-water clusters contain an average of 16 water molecules, and bulk water clusters have an average of 18 water molecules, which are adequate to approximately resemble the bulk-like environment around the central water molecule. Minor variations in the number of water molecules in a cluster do not alter the conclusions.
We sample methane-water clusters based on the distance between the methane and the nearest water molecule i.e. the central water molecule. To ensure that these clusters sample a wide variety of relevant configurations, we select the clusters with within 2.80 Å – 4.00 Å. For bulk water clusters, we tune the distance between the central water molecule and its nearest water molecule and the selected clusters have within 2.40 Å – 3.10 Å. We perform harmonic normal mode analysis based on density functional theory (DFT) calculations by using B3LYP Vosko1980 ; Lee1988 ; Becke1993 ; Stephens2002 hybrid functional in the Gaussian 16 software Frisch2016 , and 6-311++G(d,p) basis set without any geometrical optimization or any other constraints. We calculate the normal modes by constructing the Hessian matrix from the potential energy obtained by single-point energy calculations. From the eigenvalues of this Hessian matrix, we obtain the harmonic frequencies. The displacements of the atoms involved in the i-th normal mode are obtained from the eigenvectors, . As the instantaneous snapshots from the classical MD are not geometrically optimized, we obtain 7 % low-frequency range imaginary modes and are excluded from our analysis. Finally, we separate all possible real normal modes and corresponding intensities. The intensity of any normal mode is calculated using the following expression
(1) |
where represents dipole moment and is the displacement of the i-th normal mode. We evaluate the density of states (DOS), and intensity-weighted density of states (IDOS) i.e. spectral density, , by using the following equations.
(2) |
(3) |
We obtain the density of states, , by constructing a normalized histogram of vibrational frequencies of all normal modes binned to windows of 20 cm-1 and plot as a function of the bin center frequencies. We construct the spectral density from the same histogram by weighting each normal mode with the corresponding intensity.
III Results and Discussion
III.1 Full dimensional response

The instantaneous normal mode analysis provides microscopic information about the mixing of different modes, the relative symmetry of normal modes, and the extent of delocalization of normal mode vibrations Cho1994 . We represent the DOS, and the IDOS, of bulk water and methane-water clusters in Figure 2. These data include all modes, even modes from cluster surface vibrations. It is clear from Figure 2, that the DOS of bulk water reproduced almost all the features as reported earlier Biswas2017 . The normal mode analysis of only water clusters predicts the shifts in peak position and line-width relative to the experiment Biswas2017 . Despite these limited dissimilarities, one can easily recognize the qualitative similarities with the experimental spectrum. The center of the stretch region is arising at a lower frequency (3340 vs. 3400 cm-1) and wide as compared to the experimental spectrum (314 vs. 265 cm-1). The peak at 1650 cm-1 is the closest to the experimental data and is found to be subjugated by the water bending character. The librational motions are observed at a significantly higher frequency, 715 cm-1, than that observed in the experiment. These observations corroborate that DFT calculations considerably overestimate the strength of the hydrogen bonding interaction in water clusters Biswas2017 . Another difference we observe in the present case is that a bimodal distribution arises in the case of water stretching modes. This might originate from the slight difference in nuclear configurations sampled by the rigid water model as compared to the flexible water model Biswas2017 . The symmetric stretching modes show a peak maximum at 3300 cm-1, and the anti-symmetric stretching modes show a peak maximum at 3430 cm-1. In the methane-water system, we find three new signatures in the DOS. One peak arises at 1350 cm-1 because of the umbrella motion of 3H atoms about the other C-H bond. The peak around 1560 cm-1 is originating from the bend motion of the methane molecule, and another high-frequency signature arises in the frequency range of 3000–3200 cm-1, from the anti-symmetric C-H stretching mode. Besides, for the methane-water system, a decrease in intensity is also observed throughout the water stretch region compared to the bulk water data. The lower intensity is reflecting the reduction in the number of average water molecules in the case of methane-water clusters. The IDOS, of bulk water and methane-water systems are shown in Figure 2b. As reported earlier Biswas2017 , we observe a similar trend in bulk water response except for the bimodal distribution in the stretch region. In the case of methane-water response, we indeed observe different features as predicted by the DOS data. However, the respective methane-active normal mode’s low intensity gives rise to a less prominent signature in the spectral density.

Furthermore, we investigate the difference in DOS, , and that in IDOS, to probe the effect of the non-polar methane moiety on the distributions of DOS and IDOS more clearly. We obtain and by subtracting the water response from the methane-water response. In both cases, we observe negative values which signify the loss of a few water molecules in the methane-water clusters. We observe clear signatures for methane active normal modes in the case of .
In order to isolate the vibrations of the central water molecule from the contributions of surface water, we have selectively calculated the DOS and IDOS for those normal modes in which only the central water molecule has non-zero contribution. To achieve this, we track the normal mode displacement vector and isolate those normal modes in which the central water molecule has a non-zero displacement (typically 0.01 Å) (see supplementary material).
We also isolate normal modes based on the participation of the methane molecule and the central water molecule. We track the normal mode displacement of these two molecules to isolate the modes in which only the methane molecule is active and only the central water molecule is active. Only water active modes are isolated based on the displacement of the atoms in the central water molecule 0.01 Å and the displacement of the atoms in the methane molecule 0.01 Å. The only methane active modes are obtained by using the reversed condition. In Figure 3, we represent the isolated DOS and the spectral density of methane-water and bulk water systems with various conditions.

We find the corresponding DOS and IDOS of the bulk water system, and the only water active normal modes in the methane-water system are almost similar (see supplementary material). We observe three additional features in the DOS of the methane-water system, as explained earlier. By separating the only methane active modes, we notice the disappearance of the bimodal stretching frequency, which is present in bulk water, and the only water active data of methane-water clusters. Hence, it is evident that the bimodal distribution of DOS and IDOS in the high-frequency range (3100 - 3550 cm-1) originates from the water active modes. Despite an overall reduction of water active modes in the only methane active modes, it is clear from Figure 3 that there exists a small coupling of these methane active modes with the motion of other water molecules present in the cluster. To quantify the contribution of water motions in the only methane active data, we have estimated the DOS and IDOS with the displacement of the atoms in the methane molecule Å, and displacement of atoms of all other water molecules Å.

III.2 Stretch and bend modes
To investigate the mixing of stretch and bend vibrations, we decompose the spectral density and DOS by tracking the atomic stretch and bend displacements of the central water molecule for both systems. We consider a normal mode to have a bend character if the change in the H–O–H angle of the central water is . We consider a normal mode to have a stretch character if any O–H bond of the central water molecule has a displacement 0.01 Å. By comparing these two displacements, we subdivide the normal modes into a stretch only, bend only or both stretch and bend characters.
We indeed check for different cut-off values for these two displacements. Of course, the choice of the cut-offs influences the conclusions regarding the extent of mode mixing. However, several qualitative deductions remain unaffected. We select these cut-offs in such a way so that we can observe the mode mixing with higher resolution. The dependence of the degree of mode mixing on the choice of different cut-off values is represented in the supplementary material.
The decomposition of modes with only stretch and only bend characters are represented in Figure 4. This data is further decomposed based on the relative participation of the methane and the central water molecule. As explained earlier, we decompose only water active modes based on the displacement of the atoms in the central water molecule 0.01 Å and the displacement of atoms in the methane molecule 0.01 Å. The only methane active modes are obtained by using the reversed condition. In the case of both active modes, we use displacement of the atoms in the central water molecule and the methane molecule to be 0.01 Å.

It is apparent that the density of only stretch active states in both methane and water active conditions significantly changes. However, the density of only bend active states is not influenced by the presence of methane to a great extent. Interestingly, in the case of both active modes of the stretch-only distribution, the asymmetric stretch contribution is significantly suppressed, and the bimodal feature in the high-frequency range is no longer extant. This indicates that the methane active modes are coupled with symmetric stretch modes of water. To identify the nature of the methane modes, we further perform the symmetry analysis of the different normal modes, and this is presented in section C. Although the ‘bend only’ spectral density of only water active and both methane and water active modes suggests that water bend mode is not getting affected by the methane modes significantly, however, the only methane active modes have small contribution from the other water molecules as well. This suggests that the water bend mode is also coupled with the methane active modes. We also examine the extent of coupling between stretch and bend character by tracking both active modes. This is shown in Figure 4. It is evident from the figure that the stretch and bend character are indeed coupled.
III.3 Vibrational modes and symmetry
In this section, we analyze the DOS and IDOS of all the clusters for both systems with respect to the symmetry of vibration of the central water molecule. We track down the relative phase of the motion of the two O-H bonds of the central water molecule. If these O-H bonds are vibrating in phase, we call that mode symmetric mode and anti-symmetric mode otherwise. The definition of symmetric and anti-symmetric modes is shown pictorially in Figure 5. We present the results in Figure 5. The peak at 3330 cm-1 corresponds to the symmetric O-H stretching, and that at 3430 cm-1 corresponds to the asymmetric O-H stretching (see supplementary information for the bulk water data). It is also clear from the figure that both the symmetric and the asymmetric stretching modes contribute to the water bend and libration modes. However, the contribution from the asymmetric modes is a bit less than that of the symmetric modes. It is also important to realize that in both of these decomposed distributions of DOS and IDOS, we observe the signatures of the methane vibrational motions. Both of these modes also have contributions from the bend and umbrella modes of methane molecule almost equally. However, the asymmetric stretching vibration of methane is more pronounced in the symmetric stretching ensemble of the central water. To further quantify this, we investigate both active mode of methane-water data in terms of the symmetry of vibrations of the central water molecule and found that asymmetric modes do not contribute to the DOS and IDOS shown in Figure 6(see supplementary information). Only methane active data (shown in Figure 3) also suggest the same.

III.4 Delocalization of normal modes
To analyze the degree of delocalization of different normal modes, we calculate participation ratio (PR)Sastry2001 using the following equation.
(4) |
where N is the total number of atoms in the cluster, and refers to the 3N Cartesian coordinates describing the atomic displacements for each normal mode. Thus, for a local mode, the PR is of the order of , and approaches unity if it is fully delocalized over the cluster. We also calculate the degree of delocalization by estimating the number of atoms involved in each mode. The following equation calculates the number of atoms involved in the mode.
(5) |
Here is the displacement vector of the j-th atom in a given normal mode, is the distance cut-off, and H is the Heaviside step function. We calculated the participation ratio and the degree of delocalization, in units of the number of atoms, for those modes which have vibrational amplitude more than 0.01 Å using the equations 4 and 5. Figure 7 represents the PR and the degree of delocalization of different modes present in the methane-water clusters. The general observation from the figure is that the low-frequency modes are found to be more delocalized, and the high-frequency modes are localized Biswas2017 . The low frequency delocalized characters originate from the collective vibrations of the oxygen atoms in the hydrogen-bonded network. In addition to the different water vibrational modes, the figure also depicts the nature of the normal modes in which methane atoms participate actively. It is evident from the figure that the methane umbrella mode, bend mode, and asymmetric stretch modes are highly localized in nature.
IV Conclusions
It has been shown earlier that despite certain shortcomings, the cluster-based instantaneous harmonic normal mode analysis predicts different features of the experimental observations qualitatively Biswas2017 . This method is beneficial for understanding mode-mixing and interpreting the solution-phase IR spectra. In this work, we perform a detailed analysis of the methane-water system and compare the observations with bulk water response. In general, the methane-water system shows a decrease in intensity throughout the water stretch region compared to the bulk water data. This lower intensity reflects the reduction in the number of average water molecules in the case of methane-water clusters. The stretch-bend decomposition of the full DOS and IDOS reflects that only stretch active states in the presence of both methane and water active condition significantly changes, but the presence of methane does not influence the density of only bend active states to a great extent. This indicates a pronounced coupling of the methane asymmetric stretch mode with the water symmetric stretch mode. Also, we decompose the DOS and IDOS based on the symmetry of the vibrations and find that the methane active modes are coupled with symmetric stretch modes of water. We examine both active modes of methane-water data regarding the symmetry of vibrations of the central water molecule and find that asymmetric modes do not contribute to the DOS and IDOS. We investigate the extent of delocalization of different vibrational modes by estimating the participation ratio and also by identifying the number of atoms involved in any mode. This reveals that the vibrational modes having non-zero contributions from methane molecule are extremely localized in nature.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgments
RB acknowledges IIT Tirupati for support through the new faculty seed grant (NFSG), and also the computational support.
Appendix A: Supplementary Material
See supplementary information for the density of states, D() and the spectral density of states, of bulk water and methane-water with the central water molecule displacement 0.01 Å. The supplementary information also contains D() and of both bulk water and only water active case of methane-water, and their differences, and , measured as methane-water relative to bulk water. We present D() and of bulk water and methane-water systems, including the imaginary frequencies, and with different distance and angle cut-offs. Finally, We present the symmetric and anti-symmetric decomposition of modes for bulk water system.
References
- (1) C. Tanford, The Hydrophobic Effect: Formation of Micelles and Biological Membranes, Wiley-Interscience, New York, 1973.
- (2) G. Srinivas and B. Bagchi, J. Chem. Phys. 116, 8579 (2002).
- (3) G. Srinivas and B. Bagchi, J. Phys. Chem. B 107, 11768 (2003).
- (4) P. E. Mason, G. W. Neilson, J. E. Enderby, M. L. Saboungi, C. E. Dempsey, A. D. MacKerell, and J. W. Brady, J. Am. Chem. Soc. 126, 11462 (2004).
- (5) D. Chandler, Nature 437, 640 (2005).
- (6) H. J. Dyson, P. E. Wright, and H. A. Scheraga, Proc. Natl. Acad. Sci. U. S. A. 103, 13057 (2006).
- (7) A. A. Bakulin, M. S. Pshenichnikov, H. J. Bakker, and C. Petersen, J. Phys. Chem. A 115, 1821 (2011).
- (8) M. K. Hazra, S. Roy, and B. Bagchi, J. Chem. Phys. 141, 18C501 (2014).
- (9) J. A. Long, B. M. Rankin, and D. Ben-Amotz, J. Am. Chem. Soc. 137, 10809 (2015).
- (10) J. Grdadolnik, F. Merzel, and F. Avbelj, Proc. Natl. Acad. Sci. U. S. A. 114, 322 (2017).
- (11) H. S. Frank and M. W. Evans, J. Chem. Phys. 13, 507 (1945).
- (12) W. Kauzmann, Adv. Protein Chem. 14, 1 (1959).
- (13) P. J. Rossky and D. A. Zichi, Faraday Symp. Chem. Soc. 17, 69 (1982).
- (14) T. Lazaridis and M. E. Paulaitis, J. Phys. Chem. 96, 3847 (1992).
- (15) H. S. Ashbaugh and M. E. Paulaitis, J. Phys. Chem. 100, 1900 (1996).
- (16) K. A. Sharp and B. Madan, J. Phys. Chem. B 101, 4343 (1997).
- (17) G. Hummer, S. Garde, A. E. García, M. E. Paulaitis, and L. R. Pratt, J. Phys. Chem. B 102, 10469 (1998).
- (18) D. Paschek., Zeitschrift für Physikalische Chemie 223, 1091 (2009).
- (19) M. Montagna, F. Sterpone, and L. Guidoni, J. Phys. Chem. B 116, 11695 (2012).
- (20) N. Galamba, J. Phys. Chem. B 118, 4169 (2014).
- (21) B. P. H. K. D. JONG, J. E. WILSON, G. W. NEILSON, and A. D. BUCKINGHAM, Mol. Phys. 91, 99 (1997).
- (22) P. Buchanan, N. Aldiwan, A. Soper, J. Creek, and C. Koh, Chem. Phys. Lett 415, 89 (2005).
- (23) D. Laage, G. Stirnemann, and J. T. Hynes, J. Phys. Chem. B 113, 2428 (2009).
- (24) P. N. Perera, K. R. Fega, C. Lawrence, E. J. Sundstrom, J. Tomlinson-Phillips, and D. Ben-Amotz, Proc. Natl. Acad. Sci. U. S. A. 106, 12230 (2009).
- (25) K.-J. Tielrooij, J. Hunger, R. Buchner, M. Bonn, and H. Bakker, J. Am. Chem. Soc. 132, 15671 (2010).
- (26) S. Strazdaite, J. Versluis, E. H. G. Backus, and H. J. Bakker, J. Chem. Phys. 140, 054711 (2014).
- (27) B. M. Rankin, D. Ben-Amotz, S. T. van der Post, and H. J. Bakker, J. Phys. Chem. Lett. 6, 688—692 (2015).
- (28) X. Wu, W. Lu, L. M. Streacker, H. S. Ashbaugh, and D. Ben-Amotz, Angew. Chem. Int. 57, 15133 (2018).
- (29) J. Qvist and B. Halle, J. Am. Chem. Soc. 130, 10345 (2008).
- (30) Y. L. A. Rezus and H. J. Bakker, Phys. Rev. Lett. 99, 148301 (2007).
- (31) C. Petersen, K.-J. Tielrooij, and H. J. Bakker, J. Chem. Phys. 130, 214511 (2009).
- (32) L. Rossato, F. Rossetto, and P. L. Silvestrelli, J. Phys. Chem. B 116, 4552 (2012).
- (33) D. Ben-Amotz, J. Am. Chem. Soc. 141, 10569 (2019).
- (34) A. J. Bredt and D. Ben-Amotz, Phys. Chem. Chem. Phys. 22, 11724 (2020).
- (35) P. M. Wiggins, Physica A 238, 113 (1997).
- (36) K. Lum, D. Chandler, and J. D. Weeks, J. Phys. Chem. B 103, 4570 (1999).
- (37) D. M. Huang and D. Chandler, Proc. Natl. Acad. Sci. U. S. A. 97, 8324 (2000).
- (38) N. Galamba, J. Phys. Chem. B 117, 2153 (2013).
- (39) L. F. Scatena, M. G. Brown, and G. L. Richmond, Science 292, 908 (2001).
- (40) T. M. Raschke and M. Levitt, Proc. Natl. Acad. Sci. U. S. A. 102, 6777 (2005).
- (41) J. G. Davis, K. P. Gierszal, P. Wang, and D. Ben-Amotz, Nature 491, 582 (2012).
- (42) F. H. Stillinger, J. Solution Chem. 2, 141 (1973).
- (43) C. J. Fecko, J. D. Eaves, J. J. Loparo, A. Tokmakoff, and P. L. Geissler, Science 301, 1698 (2003).
- (44) K. D. Reddy and R. Biswas, J. Chem. Phys. 153, 094501 (2020).
- (45) D. Hecht, L. Tadesse, and L. Walters, J. Am. Chem. Soc. 114, 4336 (1992).
- (46) D. Hecht, L. Tadesse, and L. Walters, J. Am. Chem. Soc. 115, 3336 (1993).
- (47) B. Auer, R. Kumar, J. R. Schmidt, and J. L. Skinner, Proc. Natl. Acad. Sci. U. S. A. 104, 14215 (2007).
- (48) I. V. Stiopkin, C. Weeraman, P. A. Pieniazek, F. Y. Shalhout, J. L. Skinner, and A. V. Benderskii, Nature 474, 192 (2011).
- (49) R. Biswas, J. Furtado, and B. Bagchi, J. Chem. Phys. 139, 144906 (2013).
- (50) R. Biswas, W. Carpenter, G. A. Voth, and A. Tokmakoff, J. Chem. Phys. 145, 154504 (2016).
- (51) S. A. Corcelli, C. P. Lawrence, and J. L. Skinner, J. Chem. Phys. 120, 8107 (2004).
- (52) S. A. Corcelli and J. L. Skinner, J. Phys. Chem. A 109, 6154 (2005).
- (53) S. T. Roberts, P. B. Petersen, K. Ramasesha, A. Tokmakoff, I. S. Ufimtsev, and T. J. Martinez, Proc. Natl. Acad. Sci. U. S. A. 106, 15154 (2009).
- (54) R. Biswas, W. Carpenter, J. A. Fournier, G. A. Voth, and A. Tokmakoff, J. Chem. Phys. 146, 154507 (2017).
- (55) T. Samanta, R. Dutta, R. Biswas, and B. Bagchi, Chem. Phys. Lett. 702, 96 (2018).
- (56) M. Cho, G. R. Fleming, S. Saito, I. Ohmine, and R. M. Stratt, J. Chem. Phys. 100, 6672 (1994).
- (57) M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindah, SoftwareX 1-2, 19 (2015).
- (58) W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, J. Am. Chem. Soc. 118, 11225 (1996).
- (59) H. J. Berendsen, J. R. Grigera, and T. P. Straatsma, J. Phys. Chem. 91, 6269 (1987).
- (60) H. J. Berendsen, J. P. Postma, W. F. Van Gunsteren, A. Dinola, and J. R. Haak, J. Chem. Phys. 81, 3684 (1984).
- (61) M. Parrinello and A. Rahman, J. Appl. Phys. 52, 7182 (1981).
- (62) S. Nosé and M. L. Klein, Mol. Phys. 50, 1055 (1983).
- (63) B. Hess, H. Bekker, H. J. C. Berendsen, and J. G. E. M. Fraaije, J. Comput. Chem. 18, 1463 (1997).
- (64) T. Darden, D. York, and L. Pedersen, J. Chem. Phys. 98, 10089 (1993).
- (65) S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).
- (66) C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
- (67) A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
- (68) P. J. Stephens, F. J. Devlin, C. F. Chabalowski, and M. J. Frisch, J. Phys. Chem. 98, 11623 (2002).
- (69) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. Montgomery, J. A., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, Gaussian 16, Gaussian, Inc., Wallingford CT, USA, 2016.
- (70) S. Sastry, N. Deo, and S. Franz, Phys. Rev. E 64, 016305 (2001).
Appendix A Supplementary Material
In Figure S1, we represent the detailed analysis of the density of states, D() and spectral density of states, of bulk water and methane-water systems for all the normal modes having displacement of central water molecule 0.01 Å. We also present D() and of both bulk water and only water active case of methane-water, and their differences, and in Figure S2. The and are measured by subtracting the bulk water response from the methane-water response. In Figure S3, we present D() and of bulk water and methane-water systems, including the nearly 7% imaginary frequencies for all the clusters. Figure S4 contains the spectral density of states, for clusters with central water molecule O-H displacement cut-offs ranging between 0.01 Å - 0.05 Å and central water molecule HOH angle displacement cut-offs between 1∘ - 5∘. In Figure S5, we decompose the symmetric and anti-symmetric modes based on the central water molecule O-H displacement for bulk water clusters.




