The dynamics of water are dramatically modified upon confinement in nanoscale hydrophilic silica pores. In particular, the OH reorientation dynamics of the interfacial water are non-exponential and dramatically slowed relative to the bulk liquid. A detailed analysis of molecular dynamics simulations is carried out to elucidate the microscopic origins of this behavior. The results are analyzed in the context of the extended jump model for water that describes the reorientation as a combination of hydrogen-bond exchanges, or jumps, and rotation of intact hydrogen bonds, with the former representing the dominant contribution. Within this model, the roles of surface and dynamical heterogeneities are considered by spatially resolving the hydrogen-bond jump dynamics into individual sites on the silica pore surface. For each site the dynamics is nearly mono-exponential, indicating that dynamical heterogeneity is at most a minor influence, while the distribution of these individual site jump times is broad. The non-exponential dynamics can also not be attributed to enthalpic contributions to the barriers to hydrogen-bond exchanges. Two entropic effects related to the surface roughness are found to explain the retarded and diverse dynamics: those associated with the approach of a new hydrogen-bond acceptor and with the breaking of the initial hydrogen-bond.
I. INTRODUCTION
When water molecules are confined within nanoscale spaces, their properties and especially their dynamics become very different from those of bulk liquid water.1–4 This occurs in a wide range of situations of importance in chemistry, biochemistry, geology, and engineering. This includes, for example, heterogeneous catalysis in mesoporous materials and molecular sieves,5,6 water transport in nanofluidic devices7 as well as through porous rocks,8 water confinement within reverse micelles4,9–12 and lipid membranes,13 and permeation through membrane channels.14
One of the paradigm systems that has been considered to investigate the effects of confinement on water is the family of hydrophilic silica nanopores, which includes materials such as sol-gels,22,23 MCM-41,15–18 and Vycor glass.19–21 The dynamics of water within these silica cylindrical pores is found to be retarded with respect to bulk water, both from experiments, e.g., quasi-elastic neutron scattering (QENS)15 and optical Kerr-effect (OKE) spectroscopy,22,23 and from molecular dynamics (MD) simulations.24,25
In a prior study,24 two of us established that this dynamical perturbation affects mainly the vicinal water molecules at the pore interface, while the remaining water molecules further inside the pore display almost bulk-like properties. The hydrophilic or hydrophobic nature of the pore interface was shown to dramatically change the extent of the dynamical perturbation. However, a specific determination of the key molecular features of the silica pore that govern the dynamical perturbation is still missing. Fundamental questions remain unresolved, including, for example, the relative importances of the pore surface topology and of the different chemical groups present at the silica interface, and the enthalpic and entropic contributions to the dynamical perturbation. A further puzzle is added by the pronounced non-exponential decay in the water relaxation dynamics obtained, e.g., from OKE measurements22,23 and from simulations,24–26 which suggests a broad range of relaxation time scales covering up to two orders of magnitude.
In this paper, we combine molecular dynamics simulations and analytical modeling to analyze the reorientation and hydrogen-bond (H-bond) dynamics of water confined within amorphous silica pores. We first present a site-resolved analysis of water reorientation dynamics at the pore interface. We show that this dynamics can be fruitfully analyzed through the extended jump model, whose key features are specifically calculated next to each individual pore site. We then determine the temperature-dependence of the perturbation and its enthalpic and entropic components. We finally suggest a molecular picture which highlights the importance of the pore roughness in the magnitude of the dynamical perturbation affecting interfacial water molecules.
II. COMPUTATIONAL METHODOLOGY
A. Molecular dynamics simulations
Classical molecular dynamics simulations were carried out using the LAMMPS code.27,28 Water-filled silica pores used a rectangular periodic simulation box with dimensions of 30 Å along the pore axis and 44 Å in the perpendicular directions. The number of water molecules was determined by grand canonical Monte Carlo simulations using the Towhee program29 with a chemical potential of −30 kJ/mol, a value that gives a fully water-filled pore (441 molecules).
The SPC/E model30 was used to describe the water interactions. The parameters for the silica force field have been given previously.31,32 Intermolecular interactions were evaluated with a cutoff of 12 Å and long-range electrostatic interactions were included using three-dimensional periodic boundary conditions with an Ewald summation with a tolerance of 10−4.
Five 20 ns trajectories were propagated for water confined in the hydrophilic silica pore; the trajectories were initiated from a previously equilibrated configuration,24 and differed in the length of a further equilibration stage (which had a duration of 0, 0.5, 1.0, 1.5, or 2.0 ns for each of the five trajectories) during which velocity rescaling was used to maintain the temperature at 298 K. Each trajectory was then further equilibrated for 0.5 ns using a Nosé-Hoover thermostat33,34 with a time constant of 1 ps used to maintain the temperature before and during the data collection stage. In all cases, a 1 fs time step was used. To examine dependence on temperature, three additional 20 ns trajectories were run at both T = 285 and 315 K.
B. Analysis of hydration dynamics
To gain insight into the non-exponential dynamics, a spatially resolved analysis of hydration shell dynamics was performed, taking a similar approach as recently used to study protein hydration shells.35 The silica pore surface consists of bridging oxygen and hydroxyl (silanol and geminal) oxygen H-bond acceptor sites, hydroxyl (silanol and geminal) H-bond donor sites, and non-H-bonding silicon atoms, as illustrated in Figure 1. Note that, for the purposes of classifying site types, silanol, SiOH, and geminal, Si(OH)2, H-bonding partners were grouped together. The hydration shell was defined as containing all water OH groups that were H-bonded to a pore surface site, or whose oxygen atom lay within the first peak of the silicon-water oxygen radial distribution function (rdf). Criteria were determined from the rdfs between water oxygen or hydrogen atoms and nanopore atoms, and were as follows:
- $R_{\text{OO}^n} < 3.6$Å,$R_{\text{O}^n\text{H}} < 2.5$Å, and$\theta _{\text{HOO}^n} < 30^{\circ }$for a H-bond between a water OH group and a nanopore (bridging, silanol or geminal) oxygen H-bond acceptor site,
- $R_{\text{OO}^n} < 3.6$Å,$R_{\text{OH}^n} < 2.6$Å, and$\theta _{\text{H}^n\text{O}^n\text{O}} < 30^{\circ }$for a H-bond between water and a nanopore silanol or geminal hydrogen H-bond donor site,
RSiO < 5.1 Å for the water population in the first hydration shell not H-bonded to the pore,
where O and H denote water oxygen and hydrogen atoms and Si, On, and Hn denote nanopore silicon, oxygen, and hydrogen atoms, respectively. At each time step, each OH group in the hydration shell was assigned to a particular H-bonding surface site, or to the population of water molecules in the first hydration shell but not H-bonded to the pore. The latter population was also spatially resolved via assignment of OH groups to the nearest Si atom site. In cases of ambiguity in the assignment of an OH group, it was assigned in the priority order: acceptor site > donor site > non-H-bonded population.
The four types of sites into which the nanopore surface is decomposed: bridging oxygen H-bond acceptor sites, silanol (hydroxyl) H-bond acceptor sites, silanol (hydroxyl) H-bond donor sites, and sites that spatially resolve the water molecules which are in the first hydration shell but not H-bonded to the pore.
The four types of sites into which the nanopore surface is decomposed: bridging oxygen H-bond acceptor sites, silanol (hydroxyl) H-bond acceptor sites, silanol (hydroxyl) H-bond donor sites, and sites that spatially resolve the water molecules which are in the first hydration shell but not H-bonded to the pore.
For the calculation of quantities related to the exchange of H-bond partners via large-amplitude jumps (namely jump times and jumps angles as defined below), the following tighter H-bond criteria were used to define stable H-bond states in the Stable States Picture:36
- $R_{\text{OO}^n} < 3.0$Å,$R_{\text{O}^n\text{H}} < 2.2$Å, and$\theta _{\text{HOO}^n} < 20^{\circ }$for a H-bond between a water OH group and a nanopore (bridging, silanol or geminal) oxygen H-bond acceptor site
ROO < 2.9 Å, ROH < 2.2 Å, and θHOO < 20° for a H-bond between two water molecules.
Using these definitions, all analysis presented in this work was spatially resolved into the sub-population of OH groups next to each surface site at the time origin. All probability distributions presented were constructed by weighting each calculated value by the OH-group population next to that site. Additionally, values characterizing bulk water dynamics were extracted from the subset of water molecules which were initially farther than 7.5 Å from the nanopore surface.
All error bars reported, unless otherwise stated, were calculated using block-averaging with five blocks (based on the five trajectories) and reported at a 95% confidence level using the Student t distribution.37
III. REORIENTATION DYNAMICS IN THE HYDRATION SHELL
We follow the reorientation dynamics of the water OH bonds through the reorientation time-correlation function (tcf)
where P2 is the second Legendre polynomial and e is the unit vector along a water OH bond. The C2(t) is a good approximation to the anisotropy decay measured in time-resolved infrared pump-probe experiments.4,38
Prior QENS15, OKE,22,23 and MD24–26 studies have shown that the average water orientational relaxation for all water molecules within a hydrophilic silica pore exhibits a non-exponential decay (a result that has also been observed for other liquids22,39). Two of us showed24 that this behavior is caused by the interfacial water molecules next to the pore wall. We therefore focus here on this interfacial population. Figure 2 shows the reorientation tcf at 298 K averaged over all water molecules initially present in the nanopore hydration shell, defined as detailed in Sec. II. The tcf is non-exponential and appears to follow a power-law on the 20-200 ps interval. While it is interesting that C2(t) appears to exhibit a power-law decay, for the purposes of this paper we will focus on the more general feature that the decay is non-exponential. In the following, we investigate the origins of this non-exponentiality through a detailed analysis of the hydration shell water dynamics by site type and, ultimately, by individual site.
Water reorientation time correlation function, C2(t) in Eq. (1), at T = 298 K, for all water molecules initially H-bonded to the pore or otherwise in the hydration shell. Also shown is a power law (t−α) fit to the long-time part (t > 20 ps) of the function, with exponent α = 1.3 This exponent is slightly larger than in our previous work24 because we focus here on only the interfacial layer.
Water reorientation time correlation function, C2(t) in Eq. (1), at T = 298 K, for all water molecules initially H-bonded to the pore or otherwise in the hydration shell. Also shown is a power law (t−α) fit to the long-time part (t > 20 ps) of the function, with exponent α = 1.3 This exponent is slightly larger than in our previous work24 because we focus here on only the interfacial layer.
By performing a spatially resolved analysis of the hydration shell dynamics, as outlined in Sec. II, we calculate reorientation tcfs at 298 K for the sub-population of interfacial water assigned to each pore site. From these site-specific tcfs we extract the static distribution of reorientation times τreor for individual pore sites, based on an exponential fit to each tcf between 2 and 10 ps. This distribution is shown in Figure 3. It is apparent that the non-monoexponentiality of the C2(t) decay for all hydration shell water molecules is due to an underlying broad distribution of reorientation times for water molecules in different environments across the nanopore surface.
Probability distribution of water reorientation times in C2(t), for all water molecules initially H-bonded to the pore or in the first hydration shell, at 298 K.
Probability distribution of water reorientation times in C2(t), for all water molecules initially H-bonded to the pore or in the first hydration shell, at 298 K.
We note that the tail of this distribution (τreor > 6 ps) is—like the long-time part of the hydration shell water C2(t) decay—consistent with a power-law, although the limited number of sites on the pore surface restricts the accuracy of a fit to the distribution's tail.
IV. EXTENDED JUMP MODEL PICTURE
A. Model
In order to explore the origin of the distribution of reorientation times which gives rise to the non-monoexponentiality of the hydration shell water C2(t), we use the Extended Jump Model, a theoretical model for water reorientation introduced by one of us.36,40,41 The model is based on the underlying mechanism for water reorientation, a combination of large-amplitude angular jumps performed by an OH group when it exchanges H-bond acceptors and an additional minor contribution, the reorientation of intact H-bonds between jumps (the “frame” reorientation). Within this model, the reorientation time can be expressed as,
where
The large-amplitude angular jumps are an activated process, i.e., they pass through a transition state, and can therefore be considered to be like a chemical reaction and thus characterized by an associated reaction rate constant,
Here, nX = 1 if the OH group is engaged in a stable H-bond with acceptor X and nX = 0 otherwise; I and F are, respectively, the initial and final H-bond acceptors in the exchange process. Absorbing boundary conditions in the final state F are used, so that subsequent jumps performed by each OH group are not considered. The jump time, τ0, is then extracted via a fit of exp (−t/τ0) to 1 − Cjump(t) in the interval 1 to 7 ps.
The jump contribution to the second-order reorientation time depends on both τ0 and the magnitude of the jump angle Δθ as36
Possible H-bond acceptors in the silica pore are water oxygen atoms along with nanopore bridging and hydroxyl (silanol and geminal) oxygen atoms. For any given initial acceptor, the total jump rate constant is a sum of the rate constants for jumps to all possible final acceptors. We divide final acceptors into two groups: water oxygen atoms and pore oxygen atoms, and extract the jump time for jumps to water oxygen atoms (
and
where
B. Jump times
Using this model, we calculated the jump tcf, 1 − Cjump(t), for all water in the hydration shell, as shown in Figure 4. Like the OH reorientation tcf, C2(t), the jump tcf exhibits a non-exponential decay that occurs on a much longer time scale than for bulk water.36 While it is clear from Figure 4(b) that the jump tcf to any acceptor is not well-described by a power law, the jump tcf to water oxygen acceptors displays a power-law decay similar to that of C2(t). The difference between the two jump tcfs is due to jumps to pore acceptors that lead to a faster decay of the jump tcf but which only induce a limited reorientation of the water molecule and do not contribute significantly to the reorientation tcf. The non-exponentiality of the decays in the reorientation tcf C2(t) and in the jump tcf towards water acceptors 1 − Cjump(t) both arise from the interactions of the water molecules with the pore surface and the heterogeneity of the latter.
(a) Jump time correlation function, 1 − Cjump(t) in Eq. (3), for all water molecules initially H-bonded to the pore or H-bonded to water in the first hydration shell and jumping respectively to any acceptor (red line) and specifically to a water oxygen acceptor (blue line). (b) The long-time part of 1 − Cjump(t) on a log-log scale.
(a) Jump time correlation function, 1 − Cjump(t) in Eq. (3), for all water molecules initially H-bonded to the pore or H-bonded to water in the first hydration shell and jumping respectively to any acceptor (red line) and specifically to a water oxygen acceptor (blue line). (b) The long-time part of 1 − Cjump(t) on a log-log scale.
The pronounced non-exponential decay of the jump tcf suggests the presence of a broad distribution of interfacial relaxation dynamics. This could arise from a static disorder due to the great variety of local surface structures and to the different chemical groups exposed at the silica interface, or this could also be due to dynamical disorder as suggested in glass-forming systems42 and in bulk43 and confined44,45 supercooled water, where the relaxation dynamics next to a given site fluctuates in time.
To examine these different possibilities, we first calculate the water OH jump tcf next to each individual pore surface site in order to determine the static distribution of water H-bond exchange times at the silica interface. For each of these tcfs, the decay can now be well approximated by a single exponential. While a stretched-exponential fit exp [−(t/τ)β] of the average interfacial jump tcf (Fig. 4) yields a β exponent of 0.51 revealing a pronounced heterogeneity, the distribution of β exponents for the individual jump tcfs, shown in Figure 5, peaks above β = 0.9. This demonstrates that these tcfs are much better described by single exponentials and that the spatial decomposition has resolved most of the heterogeneity. The non-exponential relaxation of the average interface tcf is thus mostly due to a broad static distribution of jump times. Another complementary analysis of the non-exponential character of the average interface tcf considers the ratio of τ0 extracted via an exponential fit over the entire range of 1 − Cjump(t) and τ0 extracted from the integration of 1 − Cjump(t). This ratio lies between 0.9 and 1.0 for the majority of sites, and is greater than 0.8 for 97% of sites, demonstrating again that the jump tcf functions are quasi-monoexponential, and that the heterogeneity observed in the jump dynamics is principally static.
Distribution of the exponents, β, from stretched-exponential fits exp [−(t/τ)β] of the jump tcfs next to each individual pore site. The dashed line shows the β value obtained from a fit of the average jump tcf to all acceptors in Figure 4.
Distribution of the exponents, β, from stretched-exponential fits exp [−(t/τ)β] of the jump tcfs next to each individual pore site. The dashed line shows the β value obtained from a fit of the average jump tcf to all acceptors in Figure 4.
We now examine the probability distribution of these jump times, shown in Figure 6. Distributions are presented for jumps to any acceptors, only jumps to new water acceptors, and only jumps to new pore oxygen atom acceptors. The distribution of jump times for all final acceptors rises at short times, peaking at ∼5.5 ps before exhibiting a non-exponential decay at long times. Because the majority of jumps, 59%, are to final water acceptors, the distribution for only jumps to water differs only slightly from that to all acceptors. In contrast, the distribution for jumps to pore oxygen acceptors is shifted to longer times, peaking around 30 ps, but the non-exponential decay at longer jump times is not distinguishably different from that for final water acceptors within the statistics of the data. Clearly, the quantitative details of the distribution of jump times depend on the nature of the final acceptor while the qualitative features do not.
Probability distribution of jump times for water molecules initially H-bonded to the pore or in the first hydration shell but H-bonded to water, for all jumps (black circles), jumps to new water acceptors (red squares), and jumps to new pore site acceptors (green diamonds), at 298 K.
Probability distribution of jump times for water molecules initially H-bonded to the pore or in the first hydration shell but H-bonded to water, for all jumps (black circles), jumps to new water acceptors (red squares), and jumps to new pore site acceptors (green diamonds), at 298 K.
Since reorientation of an OH bond mainly occurs via jumps to a final water acceptor, while jumps to a pore acceptor lead only to partial reorientation, we concentrate on jumps to water in what follows. A deeper look into the distribution of jump times to final water acceptors is provided in Figure 7, where it is decomposed according to the initial site type. The dominant component to the distribution is for waters in the hydration shell but H-bonded to another water molecule rather than a pore site. The distribution of jump times for these jumps from one water acceptor, like the total distribution to all final water acceptors, peaks at ∼6.6 ps, corresponding to a slowdown factor of approximately 2 with respect to bulk dynamics. The distributions for other initial acceptors display different behavior as they are generally broader. However, the features are obscured by the smaller number of jump events due to the limited number of pore H-bond acceptor and donor sites. Notably, the distributions for OHs initially H-bonded to a pore atom acceptor, either a bridging or silanol oxygen, extend to longer jump times, i.e., greater than 100 ps.
Probability distribution of jump times to water (
Probability distribution of jump times to water (
The distribution of jump times to water can be further explored by examination of their spatial distribution. This is shown in Figure 8 where the nanopore structure is color-coded according to the jump times for OHs assigned to a given pore surface site. The heterogeneity of the surface is reflected in the jump times. In particular, a general feature is that sites with slow jump times tend to be more buried, less accessible to water molecules, than those with faster jump times. The factors that lead to fast H-bond exchanges for one site and slow jumps at another are explored in detail in Sec. V.
Jump times to water (
Jump times to water (
C. Jump angles
The contribution to the OH reorientation time due to hydrogen-bond exchanges includes a component involving the jump angle, Δθ, as made explicit in Eq. (4). The jump angle is defined as the OI ⋅ ⋅ ⋅ OD ⋅ ⋅ ⋅ OF angle at the transition state for the exchange of hydrogen-bond acceptors for OD−HD from OI to OF. It is therefore the angle between the average OD−HD orientations before and after the H-bond exchange. In this work we have calculated the jump angle from the configuration immediately after a jump, rather than explicit identification of the transition state,36 since we have previously found that this approach is accurate as long as the trajectory configurations are saved frequently enough.46
It is interesting to examine these jump angles in addition to the jump time, τ0, which has been discussed in detail above, since disparate values of Δθ next to the different pore sites may add to the heterogeneity in water reorientation dynamics. We have calculated the probability distribution for jump angles for waters assigned to each site on the pore surface. Then, the average jump angle for a site, j, is given from the probability distribution, Pj(Δθ), as
In the following we consider particularly the average angles for jumps that end with a new acceptor that is a water,
The average jump angle, ⟨Δθ⟩, for interfacial water molecules making jumps to a water oxygen acceptor at 298 K. Results are shown for individual sites, separated into panels for each site type defined in Figure 1.
The average jump angle, ⟨Δθ⟩, for interfacial water molecules making jumps to a water oxygen acceptor at 298 K. Results are shown for individual sites, separated into panels for each site type defined in Figure 1.
The average jump angle
V. KEY FACTORS DETERMINING THE DISTRIBUTION OF JUMP TIMES
The distribution of jump times arises because water dynamics in the pore interface is perturbed relative to water dynamics in the bulk, and this perturbation is heterogeneous across the pore surface. Specifically, the surface is both topographically and chemically heterogeneous, either or both of which may contribute to the width of the distribution. Since the jump rate constant depends on the free energy for the H-bond exchange process, a related question is whether the deviation of the free energy from that found in bulk water arises from enthalpic or entropic factors. In the following, we quantify and discuss the possible contributions to the width of the jump time distribution.
A. Enthalpic effects: Activation energies
We first consider whether the chemical heterogeneity of the pore surface gives rise to an enthalpic factor in the perturbation of jump times relative to the bulk. To that end, we calculate the activation energy for jumps to water for each surface site. The rate constant for jumps to water (
Plots of
Activation energies for jump times,
Moreover, the mapping of the activation energy values onto the pore surface presented in Figure 11 shows that there exists no clear correlation with the pore topography, i.e., the pockets and protrusions present on the pore surface. In addition, there is no correlation between the jump time perturbation extracted directly from simulation (
Activation energies for jump times (
Activation energies for jump times (
B. Entropic excluded-volume effect for the approach of a new partner
We next examine the role that the pore surface topography plays in the interfacial water jump dynamics. The perturbative effect on water dynamics of the volume excluded by the presence of a solute or surface has been explored by one of us in previous works47,51 and rationalized using a picture involving the transition state for H-bond acceptor exchange.51 The transition state configuration is defined by the positions of the reorienting water OH group and the initial and final acceptors, and in particular by the jump angle. The H-bond acceptor exchange is slowed because the approach of a new acceptor is hindered by the presence of the solute or surface. Here, we quantify this effect for water next to each pore surface site using the excluded volume slowdown factor
where F is the fraction of jump transition state locations excluded by the nanopore surface, i.e., that overlap with the volume occupied by the nanopore atoms.51 Jump transition state locations are defined using individual jumps angles for each site, as given in Sec. IV C. The predicted jump time to water, taking into account the excluded volume effect, is then given by ρVτ0, bulk51.
In Figure 12 we show the distribution of excluded volume slowdown factors for interfacial water molecules decomposed according to their initial H-bond partner. The great majority of the interfacial OH population donates an H-bond to another water molecule. These OH groups are therefore tangent to the silica wall and experience an excluded-volume slowdown factor slightly smaller than ρV = 2, as previously found for flat extended hydrophobic surfaces for example.47 Larger values, ρV > 3, are dominated by OH groups H-bonded to bridging oxygens. These OH groups are more confined since many bridging oxygen atoms are found in pockets or dips in the pore surface.
(a) Probability distribution of excluded volume factors, ρV in Eq. (9), for interfacial waters, decomposed by site type, at 298 K. (b) Same as (a) but showing only the populations H-bonded to the nanopore.
(a) Probability distribution of excluded volume factors, ρV in Eq. (9), for interfacial waters, decomposed by site type, at 298 K. (b) Same as (a) but showing only the populations H-bonded to the nanopore.
In Figure 13 we show that there is a marked correlation between the jump dynamical perturbation (
Prediction of the perturbation of jump times relative to the bulk, for jumps to water oxygen acceptors, from the excluded volume (ρV) versus the actual perturbation extracted directly from simulation (
Prediction of the perturbation of jump times relative to the bulk, for jumps to water oxygen acceptors, from the excluded volume (ρV) versus the actual perturbation extracted directly from simulation (
The local topography of the pore surface, as quantified by the distribution of ρV values, is clearly a key factor in determining the width of the distribution of jump times. This is a purely entropic factor,51 arising directly from the roughness (dips and bumps) of the pore surface. However, the excluded volume slowdown effect alone is not sufficient to fully explain the perturbation of jump dynamics in the hydration shell: in Figure 13, the jump times extracted from simulation are in general larger than those predicted by the excluded volume slowdown factor. An additional perturbative factor must play a role, as outlined in Sec. V C.
C. Entropic effect for the elongation of the initial H-bond
We now consider an additional entropic factor that can explain the differences between the calculated and predicted jump times shown in Figure 13. We concentrate on waters in the hydration shell that are not initially H-bonded to the nanopore, where there is no enthalpic effect due to the strength of the initial pore-water H-bond to consider. The missing factor can be identified by considering the radial distribution functions (rdfs) between two water oxygen atoms, once the volume has been properly normalized to consider only the fraction of the space that is accessible to water and not blocked by the pore.48 These rdfs are shown for the ten sites with fastest jump times and the ten sites with slowest jump times in Figure 14(a). It is clear from these data that the first peak in the rdf is much more pronounced for slower sites than for the faster ones. This corresponds to a larger free energy barrier for elongation of the H-bond in these slower sites.
(a) Water O–O radial distribution functions,
(a) Water O–O radial distribution functions,
Calculations of the rdfs at 285 and 315 K show that the height of this first peak is nearly independent of temperature, as shown in Figure 14(b). It is therefore essentially an entropic factor. Thus, the large slowdown of the H-bond jumps is associated with the activation entropy for stretching the initial H-bond, which is much greater in these slow sites than in the fast sites or in the bulk. This is because in the slow sites, a pair of water molecules is constrained together by the pore topography, which strongly reduces the number of accessible configurations for elongated bonds. It is this restriction on the H-bond geometry that is responsible for the entropic nature of this effect—the enthalpy of the H-bond between the pair of water molecules is not much greater than between two molecules in the bulk—and differs from the excluded volume effect considered in Sec. V B and in prior studies.47,51 The latter involves the entropy associated with the approach of the new partner, whereas this effect involves the elongation of the initial bond.
Quantitatively, this new entropic slowdown factor is determined from the free energy difference between the equilibrium H-bond length req and the transition-state configuration where this H-bond is elongated to a typical value of r‡ ≃ 3.3 Å. The perturbation factor ρE can thus be expressed as
where
Figure 14(c) quantifies the role of this new entropic term for the 20 sites considered in Figure 14(a), and compares the predicted jump times with only the excluded volume for approach of a new acceptor (as in Figure 13) to that with both entropic factors included. In the case of the fast H-bond jump times,
VI. CONCLUDING REMARKS
Using molecular dynamics simulations and analytical modeling, we have determined the molecular origins of the non-exponential reorientational relaxation of water confined in hydrophilic nanopores. We have shown that the non-exponential dynamics is due to a broad underlying distribution of rate constants for the exchange of H-bonding partners for water molecules in the nanopore hydration shell, and that this distribution is in turn due to the heterogeneity of the pore surface, and principally to the variety of local topographies at the surface. This effect is mostly entropic, acting both on the elongation of the initial H-bond in the H-bond exchange process and on the approach of the new H-bonding partner. Enthalpic effects arising from the chemical heterogeneity of the pore surface bring only a minor contribution to the distribution. We have also demonstrated that the heterogeneity observed in the water dynamics is mainly static, arising from the spatial heterogeneity of the pore surface, and that dynamical heterogeneity in the water dynamics plays only a very minor role.
A key question will then be to determine how the molecular roughness of the pore surface, which strongly depends on experimental sample preparation procedures, can modulate the width of the distribution of H-bond exchange rate constants. This is the subject of ongoing work.
ACKNOWLEDGMENTS
W.H.T. acknowledges support for this work from the Chemical Sciences, Geosciences, and Biosciences Division, Office of Basic Energy Sciences, Office of Science, U. S. Department of Energy (Grant No. DE-FG02-05ER15708). D.L. acknowledges support from the Agence Nationale de la Recherche (Grant No. ANR-2011-BS08-010-01).