Signatures of Solvation Thermodynamics in Spectra of Intermolecular Vibrations.

This study explores the thermodynamic and vibrational properties of water in the three-dimensional environment of solvated ions and small molecules using molecular simulations. The spectrum of intermolecular vibrations in liquid solvents provides detailed information on the shape of the local potential energy surface, which in turn determines local thermodynamic properties such as the entropy. Here, we extract this information using a spatially resolved extension of the two-phase thermodynamics method to estimate hydration water entropies based on the local vibrational density of states (3D-2PT). Combined with an analysis of solute-water and water-water interaction energies, this allows us to resolve local contributions to the solvation enthalpy, entropy, and free energy. We use this approach to study effects of ions on their surrounding water hydrogen bond network, its spectrum of intermolecular vibrations, and resulting thermodynamic properties. In the three-dimensional environment of polar and nonpolar functional groups of molecular solutes, we identify distinct hydration water species and classify them by their characteristic vibrational density of states and molecular entropies. In each case, we are able to assign variations in local hydration water entropies to specific changes in the spectrum of intermolecular vibrations. This provides an important link for the thermodynamic interpretation of vibrational spectra that are accessible to far-infrared absorption and Raman spectroscopy experiments. Our analysis provides unique microscopic details regarding the hydration of hydrophobic and hydrophilic functional groups, which enable us to identify interactions and molecular degrees of freedom that determine relevant contributions to the solvation entropy and consequently the free energy.

Introduction

Solvation free energies modulate driving forces in chemical and biochemical processes.[1,2] They play a critical role for the thermal stability and folding of proteins[3] in aggregation processes[4,5] and in molecular recognition.[6] While enthalpic contributions to the solvation free energy can often be understood intuitively, i.e., based on the ability of a solute to form hydrogen bonds (HBs) with water in aqueous solution, microscopic understanding and quantitative prediction of entropic changes are not straightforward. Changes in water structure or entropy have been used frequently in qualitative interpretations of experimental observations, e.g., in the context of the Hofmeister series[7] or the hydration of hydrophobic solutes.[8] However, these interpretations are often found to be misleading or even qualitatively incorrect.[9,10] This provides demand for computational tools as described here, which resolve microscopic details of the thermodynamics of solvation and connect them to experimental observables. In the following, we briefly discuss several topics related to solvation, which are or have been controversial until recently and to which the present study is able to provide novel insights.

Ions with distinct charge densities can be characterized as structure-breakers and structure-makers, which aims to describe their effect on the hydrogen bond network structure of water.[11] However, the definition and quantification of water structure in this context varies between different experiments and simulation studies. Enhanced structure can refer to properties of radial distribution functions, hydrogen bond strengths, tetrahedral order of the water hydrogen bond network, or even slow down of dynamics.[9] Likewise, ion-specific effects can be understood in terms of changes in local water entropy in the hydration shell. Here, local microscopic information can provide the key for understanding the relationships between hydration water structure, dynamics, and thermodynamics.

The unfavorable solvation free energy of hydrophobic solutes is often discussed in the context of water structure in the hydration shell. The solvation of idealized small hydrophobic solutes in water results in a decrease in the entropy due to the formation of a cavity in the solvent. This entropic effect dominates the positive solvation free energy and scales with the solute volume.[12] However, a popular interpretation of the negative solvation entropy is “iceberg”[8] or clathrate formation of water in the vicinity of hydrophobic solutes. The “iceberg” picture implies a more ice-like structure of the HB network and enhanced water–water interactions. Indeed, the HB network of water preserves its connectivity, while forming cavities that host small hydrophobic solutes (spherical radius less than ≈5 Å). An enthalpic stabilization of HBs between nearest neighbors has been observed in this case.[13] A stabilization of the water HB network in the vicinity of hydrophobic solutes may also be inferred by observations of slowed down orientational dynamics[14] and HB lifetimes.

However, a recent simulation study shows that average water–water HB interaction energies (instead of HBs between nearest neighbors) in the hydration shell of hydrophobic solutes are in fact weakly destabilized.[15] Further, enhanced lifetimes of water–water HBs in the hydration shell of hydrophobic solutes are the result of a kinetic stabilization. An accurate description is provided by an excluded volume effect, which reduces the available volume for the formation of HB rearrangement transition states in the presence of the hydrophobic solute.[16] Both results indicate that small hydrophobic solutes do not induce an enthalpic stabilization of the water HB network in their hydration shell, which would be associated with “ice-like” structures.

Additional confusion is often added when changes in the solvent structure or entropy are discussed in terms of actual thermodynamic driving forces, i.e., their contributions to the solvation free energy of hydrophobic solutes in water. Changes in the solvent structure, i.e., changes in the solvent–solvent interaction energies and corresponding effects on the entropy, have a net zero effect on the solvation free energy (or thermodynamic driving force of any other process in any solvent). This is the consequence of an exact compensation between both terms for which Ben-Naim derived an insightful and general mathematical proof, which he applied to describe solvation processes and the driving forces of hydrophobic interactions.[17] In a later study, Yu and Karplus arrived at essentially the same conclusions through a more specific investigation of a theoretical solvation process.[18] Both compensating terms contribute to calorimetrically observed total changes in enthalpy and entropy, but due to their exact cancellation, the solvation free energy is purely determined by energetic and entropic terms related to direct solute–solvent interactions. Hence, even if “iceberg” formation would occur in the hydration shell of hydrophobic solutes, which describes a change in the water structure, this would not contribute to the hydrophobic effect, i.e., an unfavorable solvation free energy.[19] Instead, the entropic cavity formation term, discussed above as the dominating contribution to the positive solvation free energy of small hydrophobic solutes, describes the entropic consequence of solute–water repulsions and is therefore not subject to such cancelations.

Regarding total solvation enthalpies and entropies (i.e., including compensating effects), it is not clear how well general concepts derived from idealized spherical solutes[12] can be translated to complex molecular shapes and geometries, as well as heterogeneous chemical properties of various functional groups of a molecular solute.[20] Here, explicit solvent simulations can provide microscopic insights into local solvent properties in order to improve our understanding of the thermodynamic consequences of solvating complex molecular solutes with multiple functional groups. In the present study, we employ molecular dynamics simulations for this purpose and investigate the entropy of water within the hydration shell of simple solutes, ranging from alkali and halide ions to rare gas atoms and small molecules with polar and nonpolar functional groups. We extend the two-phase thermodynamics (2PT)[21,22] approach to a spatially resolved version (3D-2PT), which extracts information on local molecular entropies from the vibrational density of states (VDoS). This is combined with an analysis of local contributions to the pairwise additive potential energy function used in the simulation. In addition, a complementary integral equation approach is employed to compute entropic information originating within the solute–cavity volume. This is compared with 3D-2PT results, which provide information on the cavity formation entropy indirectly through perturbed water properties in the first hydration shell.

3D-2PT complements existing methods to extract spatially resolved information on solvent thermodynamic properties, particularly the solvent entropy, using permutation reduction,[23] pair correlation functions,[24] intermolecular forces, and HB connectivity[25] or integral equation theory.[26,27] However, our approach provides the advantage of linking the observed local variations of solvent entropies directly to changes in the vibrational spectrum at far-infrared frequencies, which facilitates a thermodynamic interpretation of spectroscopic observations in this frequency range. Notably, a quantitative correlation between the thermodynamics of solvation of various alcohols in aqueous solutions and their far-infrared absorption spectrum has also been established in a recent experimental study.[28] The key idea, the extraction of entropy information from diffusive and vibrational dynamics of a liquid, is conceptually related to previous work that empirically relates the diffusion coefficient[29] or the fraction of imaginary frequency instantaneous normal modes[30] to liquid entropies.

We validate our simulations and analysis protocols by comparing total solvation enthalpies and entropies to experimental thermodynamic data for the studied solutes in water. We then analyze radial profiles of the solvent entropy in the hydration shell of monatomic solutes and interpret the observations in the light of changes in the spectrum of intermolecular vibrations. Our results allow us to localize structure-making and structure-breaking effects in ion hydration shells, as well as to identify their microscopic origins. By applying the method to small molecular solutes with an anisotropic solvation environment, we identify distinct hydration water species for polar and nonpolar functional groups characterized by their individual vibrational and thermodynamic properties. Further, separation of local solute–water and water–water interactions allows us to disentangle total solvation enthalpies and entropies into effectively contributing and canceling contributions for the individual hydration water species.

Methods

Simulation Details

Our simulations were carried out with the Gromacs 4.6.1 software.[31] Force field parameters for monovalent ions were obtained from Joung and Cheatham,[32] for rare gas atoms from Guillot and Guissani,[33] and for molecular solutes from the generalized Amber force field.[34] The TIP4P-Ew model was used[35] for water. The force field potentials were chosen to allow for a one-to-one comparison of solvation enthalpies and entropies with a recent grid cell theory (GCT) simulation study.[25] Further, the chosen force fields reproduce the experimental solvation free energies with high accuracy, as discussed in the Supporting Information (SI), Table S1. Additional simulations with the SPC/E model[36] were carried out for selected solutes for comparison with integral equation theory.

In our simulations, a single solute was kept fixed in the center of a cubic box of water with an edge length from 25 to 30 Å. As monatomic solutes, we simulated the halide anions fluoride, chloride, bromide, and iodide, the alkali cations sodium and potassium, and the rare gas atoms neon, argon, and xenon. In these cases, the radial symmetry of the isotropic environment allows for analysis of the spatial distribution of solvent properties in the environment as a function of a single variable, namely, the distance to the solute atom. In addition, the molecular solutes methanol, benzene, and N-methylacetamide were studied. Here, the solvation environment is anisotropic with a smaller number of symmetry operations. Hence, we analyze solvent properties in all three-dimensions, while using existing symmetry to improve statistics. In all simulations, long-range electrostatics were treated with the particle mesh Ewald algorithm on a 1.2 Å grid with a fourth-order interpolation.[37] Systems with a net total charge were neutralized by a uniformly distributed counter charge. Short-range electrostatic and Lennard-Jones forces were shifted to zero between 9.5 and 10.0 Å. Long-range dispersion corrections were applied for the energy and pressure. Neighbor lists were updated every 10 fs (equilibration) or 8 fs (production) with a 12 Å distance cutoff. In all simulations, a 1 fs step was used to propagate the system in time. The SETTLE algorithm[38] was used to constrain the intramolecular degrees of freedom of water molecules, while the solute atoms were kept entirely fixed.

All systems were equilibrated for 1 ns in the isobaric–isothermal (NPT) ensemble at a pressure of 1 bar and a temperature of 300 K using the Berendsen weak coupling thermostat and barostat[39] with a time constant of 1.0 ps. This was followed by production runs (80 ns for atomic solutes, 100 ns for molecular solutes) in the canonical (NVT) ensemble at 300 K using a Nosé-Hoover thermostat[40] with a time constant of 5.0 ps. During the production runs, coordinates and velocities were saved every 8 fs for the subsequent analysis.

Integral equation calculations were performed within the 1D-RISM/HNC (reference interaction site model/hypernetted chain closure)[41] approximation for the argon–water system, using numerical details as described earlier[42] with a modified SPC/E water model.[43] In order to derive solvation energies and entropies, numerical derivatives of the analytical excess chemical potential expression with respect to temperature were computed at constant pressure with ΔT = ± 1 K.[44] The temperature dependence of the water density was included based on experimental data.[45]

Spatially Resolved Solvent Enthalpy and Entropy

To analyze contributions to the solvation enthalpy and entropy of the studied solutes with spatial resolution, we first computed local solute–water and water–water potential energies on a three-dimensional cubic grid centered on the solute molecule with 32 × 32 × 32 volume elements (voxels) and a grid constant of 0.5 Å. The analysis volume and the size of the voxel elements are illustrated for the benzene molecule in Figure . Local potential energies were computed as averages over the stored simulation trajectories. For each water molecule wi, electrostatic and Lennard-Jones interaction terms with the solute s and all other water molecules wj≠i were evaluated if the minimum distance between any two atoms of the distinct molecules was less than 12 Å. Potential energy contributions were then assigned to a voxel in the analysis grid based on the position of the center of mass of water molecule wi. Notably, the volume of a single voxel is significantly smaller than the volume of a water molecule for the given grid parameters. Hence, the analysis results in near-continuous spatial distributions of the average solute–water and water–water interaction potential energies (Usw(r) and Uww(r), respectively) per water molecule after averaging over the simulated trajectories. In addition, the spatially resolved water number density nw(r) and three-dimensional pair correlation function g(r) = nw(r)/nwbulk are computed simultaneously.

Figure 1

Illustration of the cubic grid of three-dimensional volume elements (voxels) in the solvation environment of a central solute molecule, here benzene. The total analysis volume contains 32 × 32 × 32 = 32768 voxels with a grid constant of 0.5 Å for which solvent properties are computed separately. For clarity only two-dimensional arrays of voxels limiting the total analysis volume on three sides are shown.

The local contribution to the solvation enthalpy per water molecule in each voxel, ΔH(r), can be expressed aswhere the factor 1/2 in the second term compensates for double counting of water–water interactions upon spatial integration to obtain the total solvation enthalpyThe reference for the average water–water interaction potential in bulk water Uwwbulk is obtained from voxels with a minimum distance to the grid center (solute center of mass) between 9 and 11 Å for atomic and molecular solutes, respectively. The local potential energies ΔUsw/ww(r) are computed using simple real-space cut-offs for Lennard-Jones and electrostatic interactions as described above. Consequently, long-ranged interactions and interactions with the uniformly distributed counter-charge for charged systems are not taken into account, although they are included in the underlying simulation. In particular for charged systems, total solvation energies obtained from integration over ΔH(r) on the analysis grid therefore need to be corrected to include these interactions as described in the SI. The above expressions for the solvation enthalpy and its spatially resolved contributions are essentially equivalent to the corresponding procedures in grid inhomogeneous solvation theory (GIST)[24] and GCT.[25] The main distinction between each individual approach lies in the estimation of the local solvent entropy.

To spatially resolve solvent molecular entropies, we adapted the 2PT approach by Goddard and co-workers for a grid-based analysis, 3D-2PT. The 2PT approach has been originally developed to describe the absolute thermodynamic properties of bulk monatomic[21] liquids from the VDoS, I(ν), and has been subsequently extended to describe molecular liquids.[22] The VDoS of translational and rotational degrees of freedom (here per molecule) can be obtained separately from the Fourier transform of the corresponding velocity time auto correlation functionHere, vCOM describes the velocity of a water molecule center of mass (COM) with mass mw, ωrot the angular velocity of the molecule around the molecular axis k with corresponding moment of inertia I, kB is Boltzmann’s constant, and T the temperature.

Integration of Itrans(ν) and Irot(ν) over the frequencies ν provides the corresponding number of degrees of freedom (in case of water molecules: ∫0∞Itrans(ν)dν = ∫0∞Irot(ν)dν = 3). In a nutshell, the 2PT approach for rigid polyatomic molecules[22] describes decompositions of the translational and rotational VDoS of a liquid into two separate additive contributions, Itrans/rotHS(ν) and Itrans/rotHO(ν). These describe the VDoS of slow diffusive and fast vibrational degrees of freedom, respectively. The former are modeled in terms of a hard sphere (HS) fluid model constructed from the particle density and diffusion coefficient observed in the liquid (the zero frequency response of the VDoS). The VDoS of the HS model can be described by an analytical function, which includes structural relaxation processes also at nonzero frequencies. This analytic HS-model VDoS is then subtracted from the total VDoS of the liquid, leaving a remaining partial VDoS with no diffusive contributions at zero frequencies, which is then modeled analytically as a continuous set of harmonic oscillators (HO). In this approach, the HS model thus effectively describes diffusive motion and anharmonic effects that would not be included in a purely harmonic treatment.

The partition function and entropy of (1) the HS model for diffusive translational degrees of freedom and (2) the HO model for vibrations can both be expressed analytically. To express the entropy of diffusive rotational degrees of freedom, the hard sphere model is replaced by a rigid rotor (RR) expression.[22] This leads to a molecular entropy expression in terms of a sum of weighted integrals of the VDoS for the distinct degrees of freedomThe expression above does not represent a rigorous expression for the molecular entropy of liquids. Instead, it describes an interpolation between the thermodynamic properties of a HS fluid (high entropy) and a solid composed of purely harmonic oscillators (low entropy) based on the diffusivity observed in the liquid. Hence, the 2PT approach provides an estimator of the liquid entropy, which however has been shown to perform well for various liquids.[46]

Equation allows us to express the entropy of water as a sum of frequency-dependent contributions distributed over the vibrational spectrum. This is shown in Figure with illustrations of the corresponding intermolecular vibrations. Collective HB bending vibrations[47] below 100 cm–1 encode the largest fraction of the total information on the entropy. Weaker contributions result from the shoulder at 200 cm–1, which corresponds to HB stretch vibrations in the water network. For bulk TIP4P-Ew water, the VDoS for frequencies below 100, 200, 300, and 400 cm–1 contains 51%, 72%, 83%, and 88% of the total 2PT entropy, respectively. This is not surprising, as kBT at 300 K corresponds roughly to the energy required for thermal excitations of a quantum harmonic oscillator at ≈200 cm–1. Hence, vibrations at higher frequencies are mainly found in the ground state under equilibrium conditions and contribute less to the entropy, such as the broad band of water librations around 600 cm–1. This applies even more to intramolecular vibrations of water molecules beyond 1000 cm–1, which are not treated in our rigid simulation model. Intramolecular water bending and OH-stretch vibrations were previously determined to contribute only about 0.1% to the total entropy of liquid water at ambient conditions.[48]

Figure 2

Bulk water VDoS at far-infrared frequencies from 0–1000 cm–1 with illustrations of intermolecular vibrations, i.e., HB stretching modes at 200 cm–1 and librational modes between 350–1000 cm–1. Below 100 cm–1, collective HB bending vibrations are observed. The color code indicates the distribution of entropy information as obtained from eq for bulk water.

In our 3D-2PT approach, the 2PT formalism is applied to every voxel of the analysis grid described previously. For this purpose, water molecules found in a given voxel at time t are followed until t + τmax to compute ensemble-averaged velocity time autocorrelation functions for translational and rotational degrees of freedom. We used τmax = 1.6 ps to maintain information regarding the position of the water molecule, while allowing full decay of the autocorrelation function and a frequency resolution of 10 cm–1 in the Fourier transform of the time-symmetric function. The bulk water number density (instead of the local water number density) was used in the HS expressions of the 2PT approach to avoid fluctuations in the vicinity of the solute, which are attributed to the local excess chemical potential and not the actual packing density including solute atoms.

Applying the 2PT methodology to the spatially resolved velocity time correlation functions provides the spatially resolved local entropy per water molecule, S(r). The bulk water entropy Sbulk is obtained in analogy to the bulk water–water interaction potential as an average from voxels far from the solute and allows us to compute the total solvation entropy as

In summary, the combination of spatially resolved solute–solvent and solvent–solvent interaction energies with the 3D-2PT entropy analysis provides detailed insights into the spatial distribution of solvation enthalpy and entropy contributions, which also define the solvation free energy. Notably, the thermodynamic process described here corresponds to the transfer of a solute molecule with fixed position from an ideal gas phase into the solution as proposed by Ben-Naim[49] for comparison with experimental data, i.e., gas–liquid partition coefficients in case of volatile susbtances. We also employ a fixed solute orientation, thus neglecting changes in the rotational partition function of nonatomic solutes upon interaction with the solvent molecules, which causes a minor systematic error in these cases. The same applies to intramolecular solute degrees of freedom, which however can be safely ignored for the rigid solutes studied here. We point out that the total solvation free energy defined by this process are often negative even for hydrophobic solutes, albeit small in magnitude. We further note that our analysis provides information on solvation at constant pressure. This is due to the comparison of local water potential energies and entropies to bulk water (far from solute) in eqs and 5. These bulk water molecules describe the average properties of water after removal of the solute and re-equilibration of the simulated system at constant pressure p. In order to obtain solvation enthalpies, ΔHsolv = ΔUsolv + pΔV, formally the partial molar volume ΔV of the solute needs to be considered. However, the pΔV term is ignored by our approach as its magnitude is very small at ambient pressures (<0.1 kJ mol–1) for the solutes considered here. Therefore, we use ΔHsolv ≈ ΔUsolv in eqs and 2.

Results and Discussion

Comparison to Experimental Solvation Enthalpies and Entropies

To validate the 3D-2PT approach, we compare total solvation enthalpies (including cutoff corrections for charged solutes, see SI) and entropies obtained by integration of spatially resolved contributions in eqs and 5 with experimental data[25,50,51] for rare gas atoms, alkali cations, halide anions, and small molecules in Figure . The direct comparison to experimental data is meaningful as the chosen force field parameters reproduce the experimental solvation free energies with high accuracy in thermodynamic integration simulations as shown in the SI (Table S1). Where available, also alternative computational results from a GCT study[25] with the same force field parameters are shown. The corresponding numerical data are provided in Table S2 in the SI. Error bars were computed as standard deviations from the separate analysis of 20 ns trajectory fragments.

Figure 3

Direct comparison of solvation enthalpies (A) and entropies (B) of monatomic (gray: rare gases; red, cations; green: anions) and polyatomic (blue) solutes obtained from molecular simulations and experimental data from refs (50) (rare gases), (51) (ions), and (25) (molecules). Further, results obtained in this study from eqs and 5 (circles) are compared with previous GCT calcluations (crosses) using identical force field parameters.[25] The numerical data are provided in Table S2 of the SI.

The computed solvation enthalpies ΔHsolv are essentially determined by the underlying force field potential. We observe excellent agreement between experimental and simulated solvation enthalpies. The root mean squared deviations (RMSD) between our simulated and experimental solvation enthalpies is 5.8 kJ mol–1. The agreement of simulated and experimental solvation enthalpies and free energies shows that also the total solvation entropy of the studied solutes are described accurately by the simulation model. Correspondingly, these systems provide an excellent test for the ability of the 3D-2PT approach to describe the solvation entropy from a single equilibrium simulation. As the total solvation entropy is obtained as a sum over all local contributions, accurate reproduction of experimental data indicates also that the 3D-2PT expressions for local hydration water entropies are meaningful.

For the studied solutes −TΔSsolv spans a smaller range compared to ΔHsolv. In particular for the ionic species, the solvation enthalpies are large in magnitude and negative, rendering the corresponding solvation entropy seemingly less important. However, for the neutral solutes ΔHsolv and −TΔSsolv are of comparable magnitude. Hence, the entropy change cannot be neglected to describe the free energy of solvation. The smaller range of values for −TΔSsolv in Figure B allows a close look at the deviations from the experimental data. The ranking of the integrated 3D-2PT solvation entropies is predicted correctly within distinct classes of solute molecules, while the overall accuracy is comparable to GCT,[25] for example. Generally, the magnitude of the solvation entropy is slightly underestimated, except for the alkali cations, for which the opposite is the case. The RMSD for −TΔSsolv relative to the experimental data is 4.9 kJ mol–1, which can be interpreted as a sum of errors due to the underlying force field potential and systematic errors in the 3D-2PT analysis. Interestingly, experimental values for −TΔSsolv are even slightly better reproduced than the solvation enthalpies, despite the additional approximations made within the 3D-2PT approach.

Radial Profiles of Thermodynamic Solvent Properties

For atomic solutes with a spherically symmetric solvation environment, we analyze the three-dimensionally resolved solvent thermodynamic properties as a function of distance from the solute molecule. In Figure , we display radial profiles obtained for the local entropy per hydration water molecule (top panels) in the environment of anionic (left panels), cationic, and neutral atomic solutes (right panels). An analogous plot for radially resolved enthalpic contributions to the solvation free energy is shown in Figure S1 of the SI. In Figure , radially resolved water entropies are compared to radial pair correlation functions g(r) of water molecule centers of mass in the hydration environment (bottom panels) computed on the same three-dimensional grid as used for the thermodynamic analysis. The maxima of the g(r) indicate positions of distinct hydration shells, while the minima indicate interstitial regions with a decreased local water density. We use this information as a ruler to interpret the radial profiles of the local water entropy. The positions of the first maximum and minimum are indicated by arrows in each panel of Figure . Table lists the local deviations from the bulk water entropy ΔS for distances corresponding to the first maximum and minimum of the respective g(r), together with separate contributions from translational (ΔStr) and rotational (ΔSrot) degrees of freedom.

Figure 4

Radially resolved thermodynamic properties of water around monatomic solutes. Average molecular entropies of water (top panels) and water COM radial pair correlation functions g(r) (lower panels) for halide anions (left panels) and sodium and rare gas atoms (right panels). Arrows indicate positions of the first hydration shell and interstitial water, i.e., the first maximum and minimum of the respective g(r). Dashed lines for S(r) represent an extrapolation into the strongly repulsive part of the solute–water potential of mean force, where the corresponding g(r) has decreased to values below 0.2. We consider S(r → 0) = 0 as the limiting case, i.e., absence of water molecules.

Table 1

Local Deviations from the Bulk Water Molecular Entropy, ΔS, for Hydration Water at Distances rmax and rmin Corresponding to the First Maximum and Minimum of the Respective g(r) Shown in Figure a

	1stg(r) maximum				1stg(r) minimum
	brmax	cΔStr	cΔSrot	cΔS	brmin	cΔStr	cΔSrot	cΔS
F–	2.66	–10.2	–3.3	–13.5	3.42	–1.0	+0.5	–0.5
Cl–	3.16	–4.2	–1.2	–5.4	3.89	+0.9	+1.2	+2.1
Br–	3.30	–3.1	–0.7	–3.7	4.01	+1.3	+1.4	+2.7
I–	3.52	–1.6	+0.1	–1.5	4.23	+1.4	+1.4	+2.9
Na+	2.40	–11.3	–0.1	–11.4	3.28	–0.3	+0.4	+0.0
K+	2.80	–5.3	+0.1	–5.2	3.62	+0.7	+0.2	+0.9
Ne	3.30	–1.3	–0.1	–1.4	4.89	–0.0	–0.1	–0.2
Ar	3.60	–1.9	–0.1	–2.1	5.18	–0.2	–0.1	–0.3
Xe	3.92	–1.9	–0.1	–2.0	5.50	–0.2	–0.1	–0.3

Separate contributions from translational and rotational degrees of freedom are obtained according to eq . The bulk translational, rotational, and total entropies per molecule are Strbulk = 50.0 J mol–1 K–1, Srotbulk = 9.8 J mol–1 K–1, and Stotbulk = 59.8 J mol–1 K–1, respectively.

In Å

Molecular entropies in J mol–1 K–1; statistical errors within ± 0.1 J mol–1 K–1

In comparison to fluctuations in the g(r), the influence of the solute on local hydration water entropies are short-ranged. Especially for the ionic species, the maximum of the second hydration shell is clearly observed in the g(r), while local hydration water entropies are essentially converged to the bulk limit of ≈60 J mol–1 K–1 at equivalent distances. For all solutes, a decrease in the local hydration water entropy is observed in the first hydration shell relative to the bulk limit. The effect is strongest for small ions with a high charge density, i.e., fluoride and sodium, and weakest for the large iodide anion and the neutral rare gas atoms. The decrease in the local hydration water entropy in the first hydration shell is dominated by translational degrees of freedom (see Table ). Only for anions, in particular fluoride and chloride, the rotational entropy is also significantly reduced due to the formation of directional anion–water HBs.

For distances corresponding to interstitial water molecules between the first and second hydration shells, a local increase in the hydration water entropy is observed for the larger ions. In particular for chloride, bromide and iodide anions, a maximum in the radial hydration water entropy profile is observed with higher-than-bulk entropies. Both translational and rotational degrees of freedom contribute to this local entropy increase (Table ). Although weak, the effect is present also for the potassium cation, while essentially bulk-like entropies are observed for the neutral rare gas atoms at the corresponding distance from the solute. Our observations demonstrate that large ions interfere with the water HB network and disrupt its structure locally. We interpret our results as a structural transition between the first hydration shell, where water molecules are primarily bound to the ion, and the second hydration shell that has a bulk-like HB network structure. The corresponding undulations of the radial hydration water entropy profiles are also seen for the smaller fluoride and sodium ions. However, the close proximity to the ion charge center and the resulting attractive electrostatic interactions prevent higher-than-bulk entropies for interstitial water molecules in the environment of the small ions.

In a previous study, short-ranged electrostatic interactions were named responsible for ion-specific effects on solvation thermodynamics, in particular the solvation entropy.[52] Our observation of short-ranged ion-specific effects on the radial profile of the hydration water entropy for various monovalent ions confirms these findings, while providing a complete microscopic description of the respective hydration shells with a very high signal-to-noise ratio compared to previous simulation approaches.[53]

Local variations of the hydration water entropy reported in Figure are the result of changes in the VDoS of intermolecular vibrations. An analysis of these changes can provide information on the degrees of freedom and vibrational modes affected by solute–water interactions and modified water–water interactions. Furthermore, the analysis provides a connection between thermodynamic and spectroscopic properties at far-infrared frequencies, where the information content on the entropy is highest (see Figure ). In absorption[54,55] and Raman scattering[56,57] spectra of water in hydration shells, the underlying VDoS are convoluted by the corresponding IR or Raman cross sections. However, qualitative information can be obtained, e.g., based on frequency shifts of vibrational bands as discussed in the following.

The integrated intensity of the VDoS reported in Figure , for the first hydration shell (top panels) and interstitial water between the first and second hydration shells (bottom panels), corresponds to the total number of degrees of freedom per molecule and is therefore constant (see Figure S2 in the SI for detailed differences between the hydration and bulk water VDoS). Hence, a decrease in intensity at a given frequency is compensated by an increase at a different frequency. In the first hydration shell of the anions, a substantial decrease in the zero frequency response (diffusivity) and the intensity of the HB bending band are observed relative to bulk water. A new peak/shoulder arises at a higher frequency, which is attributed to the anion–water HB stretch vibration. For fluoride, this signal is most prominent and appears with ≈200 cm–1 at the highest frequency. For the heavier anions, the corresponding signal is observed between 100 and 150 cm–1. A very similar characteristic is observed for the cations, with additional peaks between 100 and 200 cm–1 attributed to vibrations of the cation–water interaction. In the quantum harmonic oscillator model, the number of thermally accessible states decreases with increasing frequency. Hence, the shift of intensity from below 100 cm–1 toward the HB stretch region at 200 cm–1 can be qualitatively assigned to a decrease in the local entropy, in agreement with the quantitative results reported in Figure and Table . The decrease in the diffusive zero frequency response within the first hydration shell of all solutes is partially also the result of the volume excluded by the solute. The latter reduces the effective volume that can be explored by nearby water molecules as the solute cavity acts as a repulsive barrier at short distances. Consequently, the 3D diffusivity as quantified by the zero frequency VDOS is effectively reduced. This represents the entropic cost of the cavity formation in our 3D-2PT analysis, which is a main contribution to the free energy of solvation of the rare gas atoms and other hydrophobic solutes as described in the Introduction. This is best observed in the top right panel of Figure S2, where it is shown that the VDoS of water in the first hydration shell of neon, argon, and xenon features a depressed zero frequency response relative to the bulk. A retardation of hydration water diffusion in the solvation shell of noble gases has also been observed experimentally, e.g., for xenon in water.[58] The decrease in zero-frequency intensity is compensated by enhanced intensities between 50 and 200 cm–1, while the remaining VDoS are essentially bulk-like. Consequently, translational contributions dominate the reduced entropy of hydration water in the first hydration shell of the rare gas atoms, as reported in Table . Reduced self-diffusion has been observed in the hydration water of many classes of solutes, including biomolecules[59−61] and hydrophobic model solutes,[62] but has not been discussed in terms of thermodynamics until recently.[63]

Figure 6

VDoS of water in the hydration environment of monatomic solutes in comparison to bulk water (thick black line). Halide anions (left panels) and alkali cations and rare gas atoms (right panels) as in Figure . Results are shown for water molecules in the first hydration shell (top panels) and interstitial waters between the first and second hydration shells (bottom panels), i.e., at distances corresponding to the first maximum and minimum of the respective g(r).

Notably, the dynamic diffusion signature of the entropic cavity formation term is obtained from properties of water molecules in the first hydration shell. This indicates a conceptual difference between simulation approaches that describe local thermodynamic properties, in particular entropies, based on particle properties in explicit solvent simulations[23−25,53] and analytical integral equation theory approaches,[4,27] which localize the cavity formation entropy contributions within the cavity volume itself. The exact cavity formation contribution to the solvation free energy can only be extracted from explicit free energy calculations, such as thermodynamic integration (see SI). It is not accessible directly in MD simulations of the final, fully coupled state in which the solute–solvent g(r) is equal to zero within the solute cavity due to repulsive interactions. Analytical theories on the other hand, such as 1D-RISM/HNC, provide finite values for the g(r) in the solute cavity, allowing the evaluation of spatially resolved solvation free energy contributions within the cavity. In Figure , we compare results obtained from a radially resolved 3D-2PT analysis and 1D-RISM calculations for the apolar solute argon. The shown cumulative integral ΔSsolvcumul(r) = ∫0ΔS(r′)nw(r′)4πr′2dr′ provides the most straightforward comparison between both approaches. We note that both curves are primarily shifted to each other. The 1D-RISM curve converges to the total solvation entropy within the solute cavity (r < 3 Å, compare to g(r) in Figure ). The 3D-2PT curve remains zero in this region and instead converges toward the full solvation entropy mostly within the first hydration shell from 3–5 Å, apart from an oscillation around 7 Å that is amplified by the r2 prefactor in the radial integral. Hence, the presence of the cavity affects the VDoS of solvating water molecules, yielding a qualitative signature of the cavity formation entropy that is detected by 3D-2PT. This is mainly relevant for nonpolar solutes, where the cavity formation entropy significantly contributes to the total solvation free energy.

Figure 5

Comparison between the cumulative integrals of radially resolved solvation entropy contributions as a function of the upper integration limit for the hydrophobic solute argon. Results are shown for a radially resolved 3D-2PT analysis of simulations in explicit solvent (solid line) and a 1D-RISM integral equation theory calculation (dashed line).

We note in passing that the solutes are immobilized in our simulation to describe a thermodynamic process of solvation, which facilitates direct comparison to experimental solvation free energies in terms of partitioning between an ideal gas and the solvent.[49] However, this treatment has an effect on the reduced mass involved in solute–solvent vibrations, which needs to be considered for comparison with experimental spectra.

Within the first hydration shell of anions, an additional effect of anion–water interactions is observed for water librational motions, i.e., hindered rotations, between 400 and 1000 cm–1. The directionality of anion–water HBs, in particular in case of fluoride, causes a blue-shift of libration frequencies. However, the entropy content of water librational modes is significantly lower due to the small number of thermally accessible states (see discussion of Figure ). Therefore, the thermodynamic consequences of the frequency shift of librational modes are less dramatic. This is well reflected by the smaller contribution of the rotational entropy to the reduced local hydration water entropy in the first hydration shell of the anions reported in Table . Cation–water interactions are less directional, and no significant effects on water libration frequencies are observed in the first hydration shell of cationic solutes.

Analysis of the VDoS for interstitial water molecules between the first and second hydration shells (g(r) minimum) allows for deeper insights into the local entropy increase observed for some of the ions. Here, the changes of the VDoS relative to bulk water are subtle, and we refer to the difference plots in Figure S2 for the detailed features. For the chloride, bromide, and iodide anions, the intensity at zero frequency is essentially bulk-like indicating bulk-like water diffusion. For the same anions, the HB bending peak intensity is increased by a few percent and features a high frequency tail up to 150 cm–1, which results in a prominent positive signal in the difference spectrum (Figure S2). As the integrated intensity of the VDoS is conserved, these increased intensities at low frequencies must be the result of red-shifted higher frequency vibrations, in particular water–water HB stretch vibrations (around 200 cm–1 in bulk water). The corresponding decrease in the VDoS intensity at this frequency can be observed in Figures and S2 but is partially compensated by an additional red-shift of the broad librational band between 300 and 1000 cm–1. Both red-shifts indicate a weakening of the water–water HBs in the interstitial region (lower force constants), which results in an increase in the local water entropy. Both effects contribute similarly to the observed entropy increase as indicated by the separate contributions of translational and rotational degrees of freedom in Table .

For the potassium cation, diffusion of interstitial water molecules is slightly enhanced compared to bulk water. The HB bending peak is bulk-like; however, additional intensity is observed in the VDoS at 100 cm–1, also due to a red-shift of water–water HB stretch vibrations. The librational modes are almost unaffected. As a result, the local hydration water entropy is increased relative to bulk. However, the effect is 2–3 times smaller than for the larger halide anions, i.e., bromide and iodide (Table ). For sodium and fluoride ions, a reduced diffusion and intensity of the HB bending vibration is also observed in the interstitial region. The resulting entropy loss is compensated by red-shifts of the water–water HB stretch vibrations and the librational modes, yielding a net zero change of the entropy relative to bulk for the sodium cation and a mild decrease in case of fluoride (Table ). The VDoS of interstitial water between the first and second hydration shells of rare gas atoms is identical to bulk water within the sampling accuracy. Consequently, the local water entropy is also equal to bulk water.

We note that the ion parameters[32] used in our simulation, while optimized to reproduce solvation free energies in TIP4P-Ew water,[35] were shown recently, together with other models, to fail in reproducing an experimentally observed increase in the water self-diffusion in solutions of cesium and potassium halides.[64] While qualitative trends for the water self-diffusion in distinct electrolytes are still reproduced, it is therefore expected that the decrease in the zero-frequency signal of the VDoS in the hydration shell of the ionic solutes is overpronounced. The correct reproduction of solvation entropies by the ion models (Figure ) suggests that short-comings in the description of hydration water diffusion might be compensated by softer, i.e., lower frequency vibrations at nonzero frequencies.

Anisotropic Distributions of Thermodynamic Solvent Properties

For molecular solutes with an anisotropic solvation environment, we apply our analysis of solvent thermodynamic properties with full spatial resolution. This allows for quantitative comparison of the effects of different functional groups on local thermodynamic properties of hydration water, while providing information on the relevant vibrational signatures in the VDoS that report on solute–water and water–water interactions. In Figure , variations of the hydration water entropy per molecule relative to bulk are shown for the hydration shell of small molecular solutes: methanol, benzene, and N-methylacetamide (NMA). Voxels with a minimum increase in the water number density by 30% relative to bulk water are displayed, i.e., the anisotropic equivalent of the first peak of radial pair correlation functions shown in Figure for atomic solutes. In Figure S3, the spatial distributions of corresponding enthalpic contributions to the solvation free energy are displayed.

Figure 7

Spatially resolved changes of the molecular entropy of water relative to the bulk in the environment of molecular solutes. Shown are voxels representing the first hydration shell with an average increase in the local water number density by at least 30% relative to the bulk, corresponding to a local excess chemical potential of < −0.5 kJ mol–1. Only voxels with a negative x-coordinate are shown for clarity. For benzene, additional voxels are shown for the π–HB site, which otherwise would be blocked from view. Voxels colors describe the local solvation entropy contribution per hydration water molecule, i.e., ΔS(r) = S(r) – Sbulk.

In order to analyze distinct hydration water species in the vicinity of various functional groups, we applied a k-means clustering algorithm[65] to identify k groups of voxels with comparable vibrational and thermodynamic properties in each solute environment. For this purpose, we simply computed the distance between two local VDoS data sets as the sum of squared intensity differences. The parameter k was determined as the maximum value at which separate VDoS with clearly distinct vibrational features were observed. Figure shows the VDoS for distinct hydration water species found for each solute in comparison to bulk water and the corresponding volume elements in the first hydration shell. Differences between the individual hydration water VDoS and bulk water are shown in Figure S4. The corresponding number of water molecules for each water species and the thermodynamic properties per water molecule are provided in Table .

Figure 8

VDoS of distinct hydration water species obtained from k-means clustering of all vibrational spectra in the hydration shell of each molecular solute. The corresponding location of voxels assigned to each water species in the three-dimensional solvation environment are identified by the corresponding color in the right-hand panels. Results are shown for methanol (A), benzene (B), and N-methylacetamide (C). Although simulations of the individual solutes are analyzed independently, the vibrational (and thermodynamic) properties of water molecules close to similar functional groups are almost identical, i.e., in the vicinity of hydrophobic solute groups and at HB binding sites.

Table 2

Average Number of Water Molecules nw and Molecular Thermodynamic Properties for Distinct Water Species Identified in the Hydration Shell of the Studied Molecular Solutes via k-means Clustering of the Local VDoS as Shown in Figure a

		noncanceling		canceling		observable
	nw	bUsw	b–TΔSsw	bΔUww/2	b–TΔSww/2	bΔH	b–TΔS
Methanol
nonpolar (far)	5.0	–1.3	+0.9	+0.5	–0.5	–0.8	+0.4
nonpolar (close)	3.4	–1.5	+1.6	+0.9	–0.9	–0.6	+0.8
interstitial	0.7	–8.2	+6.1	+5.9	–5.9	–2.3	+0.1
H-bonded	1.7	–22.4	+11.6	+8.4	–8.4	–14.0	+3.2
Benzene
nonpolar (far)	8.4	–1.7	+1.0	+0.6	–0.6	–1.1	+0.5
nonpolar (close)	6.4	–1.9	+2.3	+1.5	–1.5	–0.4	+0.8
aromatic	1.4	–12.6	+8.7	+7.9	–7.9	–4.7	+0.8
N-Methylacetamide (NMA)
nonpolar (far)	8.7	–2.0	+1.2	+0.6	–0.6	–1.4	+0.5
nonpolar (close)	5.7	–3.5	+3.2	+2.3	–2.3	–1.2	+0.9
H-bonded	2.1	–24.5	+13.4	+10.9	–10.9	–13.6	+2.4

Local free energy contributions per hydration water molecule can be obtained from the local Usw – TΔSsw (only noncanceling contributions) or ΔH – TΔS (including canceling contributions).

Free energy contributions per hydration water molecule in kJ mol–1; statistical errors within ± 0.1 kJ mol–1.

The clustering algorithm generally identified separate hydration water species in the vicinity of hydrophobic groups (far and close to the solute) and HB binding sites. In addition, a distinct water species is observed for water at the π–HB sites of the aromatic benzene ring, as well as an interstitial hydration water species in the vicinity of the methanol hydroxyl group at unfavorable HB angles. The thermodynamic and vibrational properties of each species are discussed in the following.

The vibrational and thermodynamic properties of water solvating the hydrophobic moieties of all three molecular solutes are comparable to what is observed for rare gas atoms in Figures and 6. The VDoS features a reduced zero frequency repsonse, i.e., diffusivity, compensated by increased intensites at HB bending frequencies, at ≈50 cm–1 for more distant water molecules and between 75 and 150 cm–1 for water molecules closer to nonpolar solute atoms. As described before, this VDoS signature resembles the entropic cost of cavity formation in the solvent to accommodate the solute. It results in a reduction of the molecular hydration water entropy by ≈1–3 J mol–1 K–1 (Figure ) and an unfavorable −TΔS-term of +0.4 to +0.9 kJ mol–1 per hydration water molecule (Table ). Despite the local hydration water entropy decrease, separately computed solute–water and water–water interaction energies clearly demonstrate the absence of “iceberg” formation in the vicinity of the hydrophobic groups. ΔUww is positive (Figure S3 and Table ), describing weakened water–water HBs due to a strained geometry in the presence of the solvated nonpolar groups. This is compensated by weakly attractive solute–water interactions Usw, mainly due to dispersion and molecular dipole–dipole interactions, if present. A recent comprehensive simulation study by Duboué-Dijon and Laage[15] analyzes the structure of water in the vicinity of n-propanol methyl groups using various order parameters. In particular, the tetrahedral order parameter, which accurately describes the geometry of the local HB network involving the nearest four HB partners, indicates a weak depletion of tetrahedrally structured HB configurations in the vicnity of methyl groups.[15] This finding agrees with our interpretation of changes in total water–water interaction energies in the present work.

Separate calculations of Usw and ΔUww allows us to decompose changes in the local water entropy into solute–solvent and solvent–solvent contributions. As described in the Introduction, changes in the solvent–solvent interaction energy are exactly compensated by a corresponding entropy term.[17,18] Hence, we can simply write ΔUww = TΔSww. Notably, ΔUww and TΔSww are not strictly local properties in the solute environment, as they result from pair interactions between multiple solvent molecules in distinct locations. The factor 1/2 in eq compensates for the explicit double counting of pair interactions in the calculation of total solvation energies. Therefore, we report ΔUww/2 in Figure S3 and Table as a local quantity, as well as TΔSww/2. However, the total change in the water–water interaction energy due to modified interactions with a molecule in a specific site is ΔUww. To obtain the local solute–water entropy term, we define the total change in the local solvent entropy as ΔS = ΔSsw + ΔSww/2 in analogy to eq for the enthalpy. The noncanceling solute–water energy and entropy terms Usw and −TΔSsw are then defined as strictly local contributions, which determine the solvation free energy.

The decomposition into canceling and noncanceling free energy contributions in Table shows, that contrary to an “iceberg” model, the water–water entropy increases (favorable, negative −TΔSww/2) in the vincity of hydrophobic sites as a consequence of strained water–water HBs (positive ΔUww). The overall decrease in the local solvent entropy is caused by the cavity formation term that determines ΔSsw for repulsive solute–water interactions, resulting in unfavorable, positive −TΔSsw and −TΔS terms.

For water molecules solvating the HB sites of the polar solutes methanol and NMA, we observe broken water–water HBs, as indicated by ΔUww/2 reported in Table . For water hydrating the methanol hydroxyl group ΔUww/2 amounts to +8.4 kJ mol–1. Considering that this describes only the local half of the total change in the water–water pair interaction energy, this value compares well to the average water–water HB binding energy of −17.5 kJ mol–1 in bulk water for the TIP4P water model,[66] i.e., −8.75 kJ mol–1 per particpating water molecule. The average solute–water interaction energy, ΔUsw, can be interpreted as the strength of a solute–water HB. With −22.4 kJ mol–1 and −24.5 kJ mol–1 for the methanol and NMA HB binding sites, this interaction is stronger than an average water–water hydrogen bond. The formation of the solute–water HB comes at an entropic cost quantified by −TΔSsw. This reduces the local contribution to the solvation free energy ΔUsw – TΔSsw to roughly 50% of the favorable and negative ΔUsw.

The net negative local contributions of the methanol and NMA HB binding sites to the observable solvation entropies, shown in Figure and reported by the resulting positive −TΔS terms in the far-right column of Table , represent the sums ΔSsw + ΔSww/2 and −T (ΔSsw + ΔSww/2), respectively. These reduced local hydration water entropies at HB sites of polar solutes are reflected by changes in the local VDoS shown in Figures and S4, which are qualitatively similar to water in the first hydration shell of anions. Cavity formation plus strong solute–water HBs result in a pronounced intensity decrease at zero frequency and in the HB bending band. Instead, a new peak arises between 175 and 200 cm–1, which corresponds to the solute–water HB stretch vibration. This is accompanied by a mild blue-shift of the librational band due to the directional restraint imposed by the solute–water HB.

Table also indicates broken water–water HBs at the π–HB sites on both sides of the benzene ring. This is suggested by the positive magnitude of ΔUww/2, which is close to the one observed for HB binding sites of polar solutes. We note that we report here average properties of water molecules in the π–HB binding site. Close inspection of Figure indicates that local hydration water entropies in this location are quite heterogeneous, which is not further considered here. With ΔUsw = +12.6 kJ mol–1, the interaction energy of the average π–HB is roughly 5 kJ mol–1 weaker than an average water–water HB. This estimate compares well with an enthalpic cost of 7 ± 2 kJ mol–1 for the conversion of a water–water HB into a π–HB, which was obtained from temperature-dependent Raman spectroscopy of benzene hydration water.[57] It is noted that the π–HB interaction is only described via electrostatic interactions between point charges in our force field model and hence is subject to significant oversimplification, despite this encouraging agreement.

The entropic loss due to the π–HB interaction, −TΔSsw, compensates roughly 70% of ΔUsw, i.e., (+8.7 kJ mol–1). Hence, it further reduces the binding free energy for this noncanonical interaction compared to solute–water HBs of methanol or NMA. Interestingly, the sum −T (ΔSsw + ΔSww/2) results in +0.8 kJ mol–1, which reports on the local contribution to the total solvation entropy. Despite the attractive solute–water interaction in this case, this value is comparable to what is found for water hydrating other hydrophobic groups. This observation can be traced back to the presence of low frequency librations, which we discuss in the next paragraph. In terms of local contributions to the solvation free energy, ΔUsw – TΔSsw, the π–HB site is weakly hydrophilic. However, in average only 1.4 water molecules occupy the two π–HB sites of the benzene molecule.

Nevertheless, these water molecules feature characteristic vibrational properties as shown in Figures and S4. As for other hydration water molecules that form solute–water HBs, we observe reduced intensities at zero frequency and for the HB bending term. An additional peak arises for the stretch vibration of the π–HB interaction, which is found between 100 and 150 cm–1 compared to 175 to 200 cm–1 for solute–water HBs due to the lower force constant. However, this is accompanied by a pronounced red-shift of the librational modes, which results in a pronounced peak at 400 cm–1. Hence, reduced translational entropies due to the cavity formation term and the formation of π–HB interactions are partially compensated by an increase in rotational entropy. The latter is a consequence of the less pronounced directionality of the π–HB interaction between a water molecule and the aromatic ring, which grants the bound water molecule increased rotational flexibility compared to molecules in the bulk water HB network.

As mentioned above, the k-means clustering algorithm also identified a hydration water species in the vicinity of the methanol hydroxyl group at suboptimal positions for solute–water HB formation. For these water molecules, water–water HBs seem partially broken based on values for ΔUww in Table that correspond to 60–70% of the correspodning values observed for the main solute–water HB binding site. Likewise, the solute–water interaction is only partially intact with interaction energies Usw of 25–30% of the interaction energy of an intact solute–water HB. The entropy loss due to the solute–water interaction, i.e., −TΔSsw, is essentially compensated by the local increase in the entropy due to broken water–water interactions. Hence, the entropy of these water molecules is essentially bulk-like. We characterize these water molecules as being in a transition state between water molecules that are fully integrated into the water HB network and water molecules forming a HB with the hydroxyl group of the methanol, i.e., a transition between the first and second hydration shells. We note, that only 0.7 water molecules are found on average in positions representing this transition state. Due to the strained geometries of solute–water and water–water HBs, the VDoS for these water molecules feature a broad shoulder of increased intensities between 100 and 175 cm–1 in Figures and S4, instead of a sharper peak at higher frequency. This is accompanied by a pronounced red-shift of librational modes, resulting in a local VDoS that closely resembles interstitial water molecules between the first and second hydration shells of large anions shown in Figures and S2. This shows that the almost bulk-like entropy for this hydration water species is caused by canceling effects on entropy contributions of translational and rotational degrees of freedom.

From the averaged thermodynamic properties in Table , we can compute the contribution of each hydration water species to the total solvation enthalpy and entropy by multiplying the corresponding ΔH or −TΔS with the average number of water molecules nw associated with each species. Summing over these contributions recovers ≈50–70% of the total solvation enthalpies and entropies obtained from integrating over the entire solvation environment as reported in Figure and Table S2. Consequently, the first hydration shell plays the lead role in determining the thermodynamics of solvation, while still a substantial fraction is determined from longer ranged effects. This is expected due to the treatment of water–water interaction energies and entropies, which cannot be unambiguously assigned to a single water molecule or location in the solute environment. However, while water–water interactions contribute to observable solvation enthalpies and entropies, their contributions cancel each other with respect to the solvation free energy. The latter can be expressed purely in terms of Usw and −TΔSsw, which can be assigned unambiguously to the position of the water molecule interacting with the solute. Multiplying the sum of both noncanceling free energy contributions with the respective number of water molecules and summing over the distinct hydration water species recovers 70–90% of the total solvation free energy obtained from the sum of ΔHsolv and −TΔSsolv in Table S2. Hence, short-ranged interactions of the solute with its first hydration shell dominate the solvation free energy, while longer-ranged interactions are only required for predictions of the experimentally observable solvation enthalpies and entropies, i.e., changes in the water–water interactions due to the presence of the solute.

Conclusions

In this work, we employ atomistic molecular dynamics simulations to determine thermodynamic properties of hydration water. In particular, we obtain information regarding the molecular entropy from the VDoS, i.e., the spectrum of intermolecular vibrational modes and self-diffusion. Detailed analysis of spatially resolved water entropies and vibrational frequencies provides microscopic insights into the degrees of freedom affected by the presence of the solute and their contributions to the solvation entropy. Vibrations in the far-infrared spectrum between 0–200 cm–1 or 0–6 THz provide the largest contribution to the entropic information, as expected based on the thermal accessibility of vibrationally excited states at low frequencies. Consequently, frequency shifts of vibrational bands in this frequency range or suppressed diffusive motion have substantial effects on the entropy. In the absence of attractive solute–water interactions, we can interpret reduced diffusivity as the dynamic signature of the entropic cost of cavity formation. The formed solute cavity becomes inaccessible to hydration water self-diffusion, effectively reducing the diffusivity of nearby water molecules. Consequently, a blue-shift of the zero frequency response in the VDoS is the only feature detected for the hydration water of rare gas atoms and hydrophobic moieties, where the entropic cost of cavity formation dominates the solvation free energy.[12] It is this feature in the VDoS, which allows the estimation of the entropic cost of cavity formation in our analysis. We corroborated our interpretation by comparison with an integral equation-based approach, which yields the cavity contribution to the solvation entropy from within the cavity, in which no water molecules are sampled in atomistic simulations.

When ions or hydrophilic groups with favorable interactions with water are solvated, our analysis shows that frequency shifts of intermolecular vibrations, in particular, HB bending and HB stretching, but also water librations between 300–1000 cm–1, can result in sizable entropy changes in addition to suppressed hydration water diffusion. Examples for the latter are found in the first hydration shell of small anions with a high charge density like fluoride, where the formation of strong anion–water HBs significantly restrains translational and rotational motion resulting in a reduced water entropy. In the case of the fluoride anion, a lower-than-bulk entropy is observed in the entire hydration environment, which can be interpreted as a structure-making effect.

For interstitial water between the first and second hydration shells of larger anions, weakened water–water hydrogen bonds result in a local increase in the hydration water entropy beyond the corresponding bulk value with equal contributions from translational and rotational degrees of freedom. This increased entropy can be interpreted as structure-breaking. The local maximum of the hydration water entropy also correlates with less structuring in the corresponding g(r) in Figure , i.e., less pronounced first maxima and minima, which are often considered as measures of the water structure to quantify stucture-making or structure-breaking effects. However, our radial analysis of the hydration water entropy clearly localizes structure-breaking effects in the interstitial region between the first and second hydration shells, i.e., in the transition region between water molecules directly interacting with a solvated ion (for example via ion–water HBs in the case of anions) and water molecules fully embedded in the water HB network.

Our analysis of the solvation of rare gas atoms and small hydrophobic moieties provides a clear thermodynamic picture of the process of solvation. We report contributions to observable total changes in the enthalpy and entropy, as well as a decomposition into canceling and noncanceling contributions to the free energy, related to changes in the water structure and solute–water interactions, respectively. The results demonstrate the absence of “iceberg” formation around hydrophobic solutes, despite reduced hydration water entropies. We observe an energetic destabilization of water–water HBs in the presence of hydrophobic solutes, which implies a distortion of the HB network geometry in the presence of a solute cavity. Overall, we find that the water HB network structure in hydration shells of hydrophobic moieties is indeed less ice-like than in bulk water.

The decrease in local hydration water entropies responsible for the negative solvation entropy of hydrophobic solutes is exclusively caused by the formation of the solute cavity. The effect can be described in its essence in terms of solute–water repulsions. The origin of reduced hydration water entropies is therefore similar to the origin of slowed-down hydration water dynamics, i.e., rearrangements in the water HB network, around hydrophobic solutes. The latter has been identified as a volume-exclusion effect due to the presence of the solute, which impedes the formation of transition states involved in HB rearrangements.[16] In this sense, the cavity formation indeed leaves a dynamic signature in the solvation environment, which is captured qualitatively in the 3D-2PT analysis via the reduced local VDoS at zero frequency for hydration water and translates into a negative contribution to the solvation entropy. This provides a comprehensive explanation for the comparibility of our analysis to alternative methods, which account for cavity formation contribution explicitly, such as 1D-RISM and thermodynamic integration.

Our analyis also provides insights into the hydration of aromatic rings. While the hydrophobic character of the C–H bonds dominates the solvation thermodynamics of benzene, our simulations identify π–HB sites that feature attractive interactions between water and the aromatic solute at the cost of breaking a water–water HB. A pronounced red-shift of librational modes for water molecules at the π–HB site indicates that this interaction is less directional and facilitates rotational motion compared to water–water HBs, hence resulting in an increase in the local rotational entropy.

The relation of hydration water intermolecular vibrations to thermodynamic properties discussed here allows a corresponding interpretation of spectroscopic data obtained in the far-infrared frequency range. For a one-to-one comparison of simulated VDoS and experimental spectra, absorption or Raman cross sections of vibrational modes and their changes in a solute hydration shell need to be taken into account. However, the vibrational signatures of ion–water vibrations[54,55,67] and changes of vibrational modes in the hydrogen bond bending region[61,68,69] have been reported experimentally and allow for a qualitative interpretation. An even quantitative correlation between far-infrared absorption changes and the solvation of alcohols has recently been demonstrated experimentally.[28]

For small molecular solutes, we identify distinct hydration water species with specific vibrational signatures and thermodynamic properties, which solvate different functional groups, i.e., HB sites or hydrophobic groups. This suggests that total solvation enthalpies, entropies, and free energies of a molecule may be constructed as a sum of solvation contributions of its functional groups, weighted by their number and exposure to the solvent. Indeed, this concept is used frequently for approximate determinations of protein solvation free energies based on atomic solvation parameters and solvent accessible surface areas (SASA).[70] However, such empirical approaches lack useful insights into local enthalpic and entropic contributions to the solvation free energy required for a mechanistic understanding of the solvation process. Further, detailed simulation studies have also shown that neighboring functional groups can significantly influence each other’s contribution to the total solvation free energy, i.e., in the case of the polar backbone and nonpolar side chain of hydrophobic amino acids.[71] Predictions based on solvation free energy increments for exposed functional groups are hence unlikely to be more accurate than other empirical predictors based on general solute and solvent properties, i.e., linear solvation energy relations (LSER).[72] Interesting nonadditive contributions to the solvation enthalpy and entropy are further expected when hydration water molecules interact directly with multiple functional groups, i.e., are complexed by a solute molecule. Such water molecules are known to contribute to the strongest protein–protein interactions,[73,74] and the accurate description of their thermodynamic contributions is therefore of particular interest. A similar situation is found in solvent-shared and solvent-separated ion pairs, where water molecules may be locked between two ions of opposite charge, resulting in significantly restrained degrees of freedom.[75] Explicit solvent simulation combined with a spatially resolved analysis of thermodynamic solvent properties and affected degrees of freedom via the vibrational density of states can provide detailed insights in such situations.