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.

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.

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.

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$
    ROOn<3.6 Å,
    $R_{\text{O}^n\text{H}} < 2.5$
    ROnH<2.5
    Å, and
    $\theta _{\text{HOO}^n} < 30^{\circ }$
    θHOOn<30
    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$
    ROOn<3.6 Å,
    $R_{\text{OH}^n} < 2.6$
    ROHn<2.6
    Å, and
    $\theta _{\text{H}^n\text{O}^n\text{O}} < 30^{\circ }$
    θHnOnO<30
    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.

FIG. 1.

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.

FIG. 1.

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.

Close modal

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$
    ROOn<3.0 Å,
    $R_{\text{O}^n\text{H}} < 2.2$
    ROnH<2.2
    Å, and
    $\theta _{\text{HOO}^n} < 20^{\circ }$
    θHOOn<20
    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 

We follow the reorientation dynamics of the water OH bonds through the reorientation time-correlation function (tcf)

(1)

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.

FIG. 2.

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.

FIG. 2.

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.

Close modal

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.

FIG. 3.

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.

FIG. 3.

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.

Close modal

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.

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,

(2)

where

$\tau ^{jump}_{reor}$
τreorjump is the contribution to OH reorientation due to the jumps between H-bond acceptors and
$\tau ^{frame}_{reor}$
τreorframe
is that due to the intact H-bond rotation. The
$\tau ^{jump}_{reor}$
τreorjump
contribution, which depends on both the time scale and magnitude of the H-bond jumps, dominates and henceforth we will only consider it.

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,

$\tau _0^{-1}$
τ01⁠. The inverse of this rate constant for H-bond acceptor exchange is the jump time, τ0, which can be obtained from the stable states cross time correlation function,

(3)

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 (−t0) 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 

(4)

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 (

$\tau _{0,\text{H}_2\text{O}}$
τ0,H2O⁠) and for jumps to pore oxygen atoms (τ0, pore) via the two relationships

(5)

and

(6)

where

$P_{\text{H}_2\text{O}}^\infty$
PH2O and
$P_{pore}^\infty$
Ppore
are the population of stable states F = O and F = On at infinite time. In practice, this corresponds to the fraction of all jumps whose final H-bond acceptor is a water oxygen or a nanopore surface site, respectively, when Cjump(t) has plateaued at Cjump(t) = 1 or can be extrapolated to 1, i.e., when all OH groups contributing to the correlation function have undergone H-bond exchange.

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.

FIG. 4.

(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.

FIG. 4.

(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.

Close modal

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.

FIG. 5.

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.

FIG. 5.

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.

Close modal

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.

FIG. 6.

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.

FIG. 6.

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.

Close modal

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.

FIG. 7.

Probability distribution of jump times to water (

$\tau _{0,\text{H}_2\text{O}}$
τ0,H2O⁠) in the nanopore hydration shell decomposed according to site type as defined in Figure 1, at 298 K.

FIG. 7.

Probability distribution of jump times to water (

$\tau _{0,\text{H}_2\text{O}}$
τ0,H2O⁠) in the nanopore hydration shell decomposed according to site type as defined in Figure 1, at 298 K.

Close modal

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.

FIG. 8.

Jump times to water (

$\tau _{0,\text{H}_2\text{O}}$
τ0,H2O⁠) in the nanopore hydration shell, mapped onto the nanopore surface, at 298 K. The two views show the two halves of the cylindrical nanopore. Also shown are the silica framework (yellow) and the hydration shell water molecules.

FIG. 8.

Jump times to water (

$\tau _{0,\text{H}_2\text{O}}$
τ0,H2O⁠) in the nanopore hydration shell, mapped onto the nanopore surface, at 298 K. The two views show the two halves of the cylindrical nanopore. Also shown are the silica framework (yellow) and the hydration shell water molecules.

Close modal

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

(7)

In the following we consider particularly the average angles for jumps that end with a new acceptor that is a water,

$\langle \Delta \theta \rangle _{j,\text{H}_2\text{O}}$
Δθj,H2O⁠. The results for all sites are shown in Figure 9, separated by initial H-bond acceptor type. We note that the jump angle distribution for molecules outside the hydration shell (not shown) is consistent with that of bulk water, with an average of 70.2°.

FIG. 9.

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.

FIG. 9.

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.

Close modal

The average jump angle

$\langle \Delta \theta \rangle _{j,\text{H}_2\text{O}}$
Δθj,H2O varies between sites, but the degree to which it does so depends on the site type. For OHs that are initially hydrogen-bonded to a pore site, the jump angle is generally larger, ranging from 66° to 95° if the initial acceptor is a bridging oxygen and 68°–99° if it is an oxygen on a silanol group. For these sites where the OH is initially pointed toward the pore surface, the jump angle to a new, water acceptor is by necessity larger. This is consistent with prior work on the water jump mechanism next to a flat extended surface47 and within the cages of a zeolite.48 The variation reflects the heterogeneity of these surface acceptor sites, particularly volume excluded by the neighboring pore atoms. In comparison, the average angles for OHs on a water accepting a hydrogen bond from a pore silanol group are between 65° and 73°. Similarly, for hydration shell waters that are not hydrogen bonded to a site on the surface, the average jump angles are, with one exception, between 61° and 72°. These values are closer to that for bulk water, indicating weaker perturbations on the geometries of hydrogen-bond exchanges.

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.

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 (

$k_{0,\text{H}_2\text{O}}\break = 1/\tau _{0,\text{H}_2\text{O}}$
k0,H2O=1/τ0,H2O⁠) and the corresponding activation energy (
$E_{a,\text{H}_2\text{O}}$
Ea,H2O
) are related by

(8)

Plots of

$\ln (k_{0,\text{H}_2\text{O},T}) [T=285, 298, 315$
ln(k0,H2O,T)[T=285,298,315 K] versus 1/T for each site are linear or approximately linear. We therefore extract the activation energy for jumps to water for each nanopore site from the slope of each Arrhenius plot. The activation energies are shown in Figure 10, and vary little by site, fluctuating around the bulk water value of 3.5 kcal mol−1.36 There clearly exists no correlation between the value of the activation energy and the site type (i.e., whether the initial H-bond acceptor is a bridging, hydroxyl, or water oxygen).

FIG. 10.

Activation energies for jump times,

$E_{a,\text{H}_2\text{O}}$
Ea,H2O in Eq. (8), for interfacial water molecules making jumps to a water oxygen acceptor. Results are shown for individual sites, separated into panels for each site type defined in Figure 1.

FIG. 10.

Activation energies for jump times,

$E_{a,\text{H}_2\text{O}}$
Ea,H2O in Eq. (8), for interfacial water molecules making jumps to a water oxygen acceptor. Results are shown for individual sites, separated into panels for each site type defined in Figure 1.

Close modal

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 (

$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk}$
τ0,H2O/τ0,bulk⁠, where τ0, bulk is the jump time in bulk water) and the perturbation which would arise from the calculated activation energies (
$\exp [-(E_{a,bulk}-E_{a,\text{H}_2\text{O}})/k_BT]$
exp[(Ea,bulkEa,H2O)/kBT]
). We therefore rule out a strong enthalpic effect in the jump time perturbation. This is unsurprising since the H-bond acceptor sites are all oxygen atoms with fairly similar partial charges (−0.64 for bridging oxygens, −0.74 for hydroxyl oxygens, and −0.8476 for water oxygens in the force fields employed for this study). This observation is in contrast to the situation in more chemically heterogenous systems, for example, protein or nucleic acid hydration shells, which contain both strong H-bond acceptors such as carboxylate or phosphate groups and weaker oxygen and nitrogen H-bond acceptors.35,49,50

FIG. 11.

Activation energies for jump times (

$E_{a,\text{H}_2\text{O}}$
Ea,H2O⁠) for interfacial water molecules making jumps to a water oxygen acceptor, mapped onto the nanopore surface. Also shown is the silica framework in yellow.

FIG. 11.

Activation energies for jump times (

$E_{a,\text{H}_2\text{O}}$
Ea,H2O⁠) for interfacial water molecules making jumps to a water oxygen acceptor, mapped onto the nanopore surface. Also shown is the silica framework in yellow.

Close modal

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

(9)

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.

FIG. 12.

(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.

FIG. 12.

(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.

Close modal

In Figure 13 we show that there is a marked correlation between the jump dynamical perturbation (

$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk}$
τ0,H2O/τ0,bulk⁠) extracted from simulation and the perturbation predicted by the excluded volume slowdown factor (ρV). This is particularly the case for jumps between two water acceptors, data for which are shown in Figures 13(c) and 13(d), which make up more than half (58%) of all jumps to water in the first hydration shell. (Recall that, as outlined in Sec. II, water OH groups assigned to hydroxyl donor sites, Figure 13(c), are initially accepting a H-bond from the pore and donating a H-bond to another water molecule.) For jumps from a pore acceptor to a water acceptor, Figures 13(a) and 13(b), the correlation between the excluded volume factor ρV and the actual slowdown
$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk}$
τ0,H2O/τ0,bulk
is much less marked, due to the additional enthalpic effect which the varying strength of the initial pore-water H-bond, relative to a water-water H-bond, has been shown to have on the jump dynamics.50 

FIG. 13.

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 (

$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk}$
τ0,H2O/τ0,bulk⁠) at 298 K. Results are separated into four panels according to site type defined in Figure 1; the black lines indicate perfect correlation. (Note τ0, bulk = 3.3 ps.)

FIG. 13.

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 (

$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk}$
τ0,H2O/τ0,bulk⁠) at 298 K. Results are separated into four panels according to site type defined in Figure 1; the black lines indicate perfect correlation. (Note τ0, bulk = 3.3 ps.)

Close modal

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.

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.

FIG. 14.

(a) Water O–O radial distribution functions,

$g_{\text{O--O}}^{pore}(r)$
gO–Opore(r)⁠, for selected fast (blue lines) and slow (red lines) individual sites for hydration shell waters not H-bonded to the pore (see the text) (b) The
$g_{\text{O--O}}^{pore}(r)$
gO–Opore(r)
for 3 selected slow sites (green, blue, and red lines) and 1 fast site (black line) at 285 K (dotted lines), 298 K (solid lines), and 315 K (dashed lines). (c) Prediction of the perturbation of jump times compared to that determined directly from the simulation (horizontal axis); predictions are shown for the entropic factor for the approach of the new H-bond partner alone (ρV) (blue) and for its product with the entropic factor for elongation of the initial H-bond (ρE) (magenta).

FIG. 14.

(a) Water O–O radial distribution functions,

$g_{\text{O--O}}^{pore}(r)$
gO–Opore(r)⁠, for selected fast (blue lines) and slow (red lines) individual sites for hydration shell waters not H-bonded to the pore (see the text) (b) The
$g_{\text{O--O}}^{pore}(r)$
gO–Opore(r)
for 3 selected slow sites (green, blue, and red lines) and 1 fast site (black line) at 285 K (dotted lines), 298 K (solid lines), and 315 K (dashed lines). (c) Prediction of the perturbation of jump times compared to that determined directly from the simulation (horizontal axis); predictions are shown for the entropic factor for the approach of the new H-bond partner alone (ρV) (blue) and for its product with the entropic factor for elongation of the initial H-bond (ρE) (magenta).

Close modal

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

(10)

where

$\Delta G^{\ddag }_{elong,bulk} = G_{bulk}(r^{\ddag }) - G_{bulk}(r_{eq})$
ΔGelong,bulk=Gbulk(r)Gbulk(req) (i.e., the barrier for the elongation of the initial H-bond in the bulk) and
$g_{\text{OO}}^{pore}$
gOOpore
and
$g_{\text{OO}}^{bulk}$
gOObulk
are the O–O rdfs next to the pore site under consideration and in the bulk, respectively.

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,

$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk} \sim 1$
τ0,H2O/τ0,bulk1⁠, the excluded volume factor alone gives a good prediction of the dynamical perturbation and there is little additional effect. However, for the sites with slower jump times,
$\tau _{0,\text{H}_2\text{O}}/\tau _{0,bulk} > 1$
τ0,H2O/τ0,bulk>1
, this entropic factor increases the predicted jump time in every case, yielding better agreement with the slowdown (relative to the bulk value) from the MD simulations. The respective RMS errors for all sites shown are 2.8 for the ρV factor only and 1.5 for the ρVρE product. The systematic increase in the predicted values for slower jump times is consistent with the missing factor evident from Figure 13, where the agreement with the calculated results are adequate for faster jump times but increasingly less so for slower jump times.

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.

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).

1.
J. C.
Rasaiah
,
S.
Garde
, and
G.
Hummer
,
Annu. Rev. Phys. Chem.
59
,
713
(
2008
).
2.
M.-C.
Bellissent-Funel
,
Eur. Phys. J. E - Soft Matter
12
,
83
(
2003
).
3.
N.
Giovambattista
,
P. J.
Rossky
, and
P. G.
Debenedetti
,
Annu. Rev. Phys. Chem.
63
,
179
(
2012
).
4.
M. D.
Fayer
and
N. E.
Levinger
,
Annu. Rev. Anal. Chem.
3
,
89
(
2010
).
5.
A.
Corma
,
S.
Iborra
, and
A.
Velty
,
Chem. Rev.
107
,
2411
(
2007
).
6.
7.
H. A.
Stone
,
A. D.
Stroock
, and
A.
Ajdari
,
Annu. Rev. Fluid Mech.
36
,
381
(
2004
).
8.
B.
Berkowitz
,
Adv. Water Res.
25
,
861
(
2002
).
9.
P. A.
Pieniazek
,
Y.-S.
Lin
,
J.
Chowdhary
,
B. M.
Ladanyi
, and
J. L.
Skinner
,
J. Phys. Chem. B
113
,
15017
(
2009
).
10.
R.
Costard
,
N. E.
Levinger
,
E. T. J.
Nibbering
, and
T.
Elsaesser
,
J. Phys. Chem. B
116
,
5752
(
2012
).
11.
A. A.
Bakulin
,
D.
Cringus
,
P. A.
Pieniazek
,
J. L.
Skinner
,
T. L. C.
Jansen
, and
M. S.
Pshenichnikov
,
J. Phys. Chem. B
117
,
15545
(
2013
).
12.
A. V.
Martinez
,
L.
Dominguez
,
E.
Małolepsza
,
A.
Moser
,
Z.
Ziegler
, and
J. E.
Straub
,
J. Phys. Chem. B
117
,
7345
(
2013
).
13.
E.
Yamamoto
,
T.
Akimoto
,
Y.
Hirano
,
M.
Yasui
, and
K.
Yasuoka
,
Phys. Rev. E
87
,
052715
(
2013
).
14.
S.
Berneche
and
B.
Roux
,
Nature
414
,
73
(
2001
).
15.
A.
Faraone
,
L.
Liu
,
C.-Y.
Mou
,
P. C.
Shih
,
J. R. D.
Copley
, and
S.-H.
Chen
,
J. Chem. Phys.
119
,
3963
(
2003
).
16.
P.
Smirnov
,
T.
Yamaguchi
,
S.
Kittaka
,
S.
Takahara
, and
Y.
Kuroda
,
J. Phys. Chem. B
104
,
5498
(
2000
).
17.
S.
Takahara
,
N.
Sumiyama
,
S.
Kittaka
,
T.
Yamaguchi
, and
M.-C.
Bellissent-Funel
,
J. Phys. Chem. B
109
,
11231
(
2005
).
18.
S.
Kittaka
,
S.
Takahara
,
H.
Matsumoto
,
Y.
Wada
,
T. J.
Satoh
, and
T.
Yamaguchi
,
J. Chem. Phys.
138
,
204714
(
2013
).
19.
F.
Bruni
,
M. A.
Ricci
, and
A. K.
Soper
,
J. Chem. Phys.
109
,
1478
(
1998
).
20.
J.-M.
Zanotti
,
M.-C.
Bellissent-Funel
, and
S.-H.
Chen
,
Phys. Rev. E
59
,
3084
(
1999
).
21.
H.
Thompson
,
A. K.
Soper
,
M. A.
Ricci
,
F.
Bruni
, and
N. T.
Skipper
,
J. Phys. Chem. B
111
,
5610
(
2007
).
22.
R. A.
Farrer
and
J. T.
Fourkas
,
Acc. Chem. Res
36
,
605
(
2003
).
23.
A.
Scodinu
and
J. T.
Fourkas
,
J. Phys. Chem. B
106
,
10292
(
2002
).
24.
D.
Laage
and
W. H.
Thompson
,
J. Chem. Phys.
136
,
044513
(
2012
).
25.
A. A.
Milischuk
and
B. M.
Ladanyi
,
J. Chem. Phys.
135
,
174709
(
2011
).
26.
A. A.
Milischuk
,
V.
Krewald
, and
B. M.
Ladanyi
,
J. Chem. Phys.
136
,
224704
(
2012
).
27.
S. J.
Plimpton
,
J. Comput. Phys.
117
,
1
(
1995
).
28.
See http://lammps.sandia.gov for the LAMMPS molecular dynamics package.
29.
See http://towhee.sourceforge.net for the Towhee Monte Carlo simulation code.
30.
H. J. C.
Berendsen
,
J. R.
Grigera
, and
T. P.
Straatsma
,
J. Phys. Chem.
91
,
6269
(
1987
).
31.
T. S.
Gulmen
and
W. H.
Thompson
,
Dynamics in Small Confining Systems VIII
, edited by
J. T.
Fourkas
,
P.
Levitz
,
R.
Overney
, and
M.
Urbakh
(
Mater. Res. Soc. Symp. Proc.
,
Warrendale, PA
,
2005
), Vol.
899E
.
32.
T. S.
Gulmen
and
W. H.
Thompson
,
Langmuir
22
,
10919
(
2006
).
34.
W. G.
Hoover
,
Phys. Rev. A
31
,
1695
(
1985
).
35.
E.
Duboué-Dijon
,
A. C.
Fogarty
, and
D.
Laage
,
J. Phys. Chem. B
118
,
1574
(
2014
).
36.
D.
Laage
and
J. T.
Hynes
,
J. Phys. Chem. B
112
,
14230
(
2008
).
37.
D. P.
Shoemaker
,
C. W.
Garland
, and
J. W.
Nibler
,
Experiments in Physical Chemistry
(
McGraw-Hill
,
New York
,
1989
).
38.
Y. S.
Lin
,
P. A.
Pieniazek
,
M.
Yang
, and
J. L.
Skinner
,
J. Chem. Phys.
132
,
174505
(
2010
).
39.
C. D.
Norton
and
W. H.
Thompson
,
J. Phys. Chem. B
118
,
8227
(
2014
).
40.
D.
Laage
and
J. T.
Hynes
,
Science
311
,
832
(
2006
).
41.
D.
Laage
,
G.
Stirnemann
,
F.
Sterpone
,
R.
Rey
, and
J. T.
Hynes
,
Annu. Rev. Phys. Chem.
62
,
395
(
2011
).
42.
R.
Zwanzig
,
Acc. Chem. Res.
23
,
148
(
1990
).
43.
G.
Stirnemann
and
D.
Laage
,
J. Chem. Phys.
137
,
031101
(
2012
).
45.
G. B.
Suffritti
,
P.
Demontis
,
J.
Gulín-González
, and
M.
Masia
,
J. Phys.: Condens. Matter
26
,
155103
(
2014
).
46.
A. A.
Vartia
,
K. R.
Mitchell-Koch
,
G.
Stirnemann
,
D.
Laage
, and
W. H.
Thompson
,
J. Phys. Chem. B
115
,
12173
(
2011
).
47.
G.
Stirnemann
,
P. J.
Rossky
,
J. T.
Hynes
, and
D.
Laage
,
Faraday Discuss.
146
,
263
(
2010
).
48.
A. C.
Fogarty
,
F.-X.
Coudert
,
A.
Boutin
, and
D.
Laage
,
ChemPhysChem
15
,
521
(
2014
).
49.
F.
Sterpone
,
G.
Stirnemann
, and
D.
Laage
,
J. Am. Chem. Soc.
134
,
4116
(
2012
).
50.
F.
Sterpone
,
G.
Stirnemann
,
J. T.
Hynes
, and
D.
Laage
,
J. Phys. Chem. B
114
,
2083
(
2010
).
51.
D.
Laage
,
G.
Stirnemann
, and
J. T.
Hynes
,
J. Phys. Chem. B
113
,
2428
(
2009
).