Illustration of the simulated system. Throughout this work, unless otherwise stated, oxygen atoms are colored red, hydrogen white, and silicon yellow. (a) Two identical silica plates were solvated in an initially 5.5 nm cubic box of SPC/E water, and immobilized at a particular separation in the x-direction to form a planar gap of fixed width. The vacant regions in this x,y-projection illustrate the position of the plates, between which resides a slab of confined water. (b) and (c) Tranverse views of the plates. For reasons described in the main text, the gap width, d, is defined as 0.066 nm greater than the linear separation between the plane containing the innermost layer of O atoms in each plate. (d) Plan view of the plates. The gap region was discretized into 42 cubic cells of dimensions [Δy = 0.495 nm, Δz = 0.429 nm], which span the gap in the x-direction. A coarse-grained characterization of a particular system configuration was developed for the purposes of linear dimensionality reduction by enumerating the water molecules residing within each cell, to synthesize what we term the occupancy vector of the gap.
P-d phase diagram for water confined between atomistically detailed hydrophobic plates at T = 298 K. Seven different phases of confined water were observed: L—liquid (cyan circles), ML—monolayer liquid (purple circles), BI—bilayer ice (red squares), MI—monolayer ice (blue squares), THF—trilayer heterogeneous fluid (pink triangles), fTHF—frustrated trilayer heterogeneous fluid (green triangles), and V—vapor (black crosses). In the vapor phase region, the + symbols correspond to an effective vacuum between the plates, whereas the × symbols correspond to state points in which cavitation was observed in the fluid surrounding the plates, indicating thermodynamic instability of the bulk liquid under these conditions. The structural and dynamical behavior of each “wet” morphology was characterized in detail at one particular state point indicated by the open squares in the diagram. Superposition of two phase-characterizing symbols (e.g., red squares (BI) and green triangles (fTHF)) indicates the observation of phase coexistence in simulations conducted at this state point. The solid lines demarcate a tentative phase map, but are intended primarily as a guide to the eye rather than a precise localization of the phase boundaries. The dashed extension of these lines indicates projected extrapolations into regions of P-d space where simulations were not conducted, or the results were too ambiguous to clearly delineate the phase behavior. Employing the contact angle of θ c = 105° computed for this system by Giovambattista et al.,102 the blue dotted-dashed line represents the vapor-liquid phase boundary predicted by Eq. (2) over the range of plate separations for which we observe a vapor-liquid phase transition. Treating the contact angle as an adjustable parameter, the predictions of Eq. (2) show relatively good agreement with the observed boundary for θ c = 130°, as illustrated by the red dotted-dashed line. The bilayer ice to vapor transition at [P ≈ 0.093 GPa, d = 0.60 nm] was investigated by the analysis of decompression-induced sublimation across this phase boundary.
(a) P-d-σ phase diagram for water confined between atomistically detailed hydrophobic plates at T = 298 K. σ, the number of water molecules per unit area residing within the gap, was computed by averaging over the final 2 ns of each 5 ns molecular dynamics simulation. The pairing of phases with symbols is identical to that in Fig. 2. For clarity of exposition in the σ dependence of the phase behavior, state points at which bulk cavitation or phase coexistence were observed are omitted in the present figure. Additionally, the state points at [P = 0.50 GPa, d = 0.45 nm] and [P = 0.60 GPa, d = 0.45 nm] were removed due to imperfections in the bilayer crystals which did not heal on the time scale of the simulations (cf. Sec. III A), resulting in artificially low σ values. The standard deviation in σ is less than ±1 nm−2 for all state points. Points sharing common d values are connected by lines as a guide to the eye. (b) Two-dimensional projections of panel (a) into the d-σ plane illustrates that for plate separations d ⩽ 0.65 nm, σ is a unique discriminant of the phase behavior. At larger plate separations, this ceases to be true. For the monolayer and bilayer ice phases, we observe good agreement with the theoretical σ values for idealized infinite crystals of 9.15 nm−2 and 18.31 nm−2. The dispersion around these values arises from edge effects in the finite sized plates employed in this work, and imperfections in the crystal structures.
P O, gap (x), for each of the six “wet” confined morphologies. In this, and the following two figures, calculations for each phase were conducted at the boxed state points in the phase diagram presented in Fig. 2, and we follow the pairing of colors and phases adopted therein: L—liquid [P = 0.20 GPa, d = 0.70 nm], ML—monolayer liquid [P = 0.20 GPa, d = 0.50 nm], BI—bilayer ice [P = 0.20 GPa, d = 0.55 nm], MI—monolayer ice [P = 0.20 GPa, d = 0.45 nm], THF—trilayer heterogeneous fluid [P = 0.60 GPa, d = 0.80 nm], and fTHF—frustrated trilayer heterogeneous fluid [P = 0.75 GPa, d = 0.65 nm]. x = 0 nm is defined as the location of y,z-mid-plane lying equidistant between the plates. Normalization of P O, gap (x) by the area of the plates, A = (3.218 × 3.215) nm2 = 10.35 nm2, yields the local density profile across the gap, ρ O, gap (x).
g O, gap (r), for each of the six “wet” confined morphologies computed at the boxed state points in the phase diagram in Fig. 2. r is the Euclidean distance between O atom centers. Normalization of g O, gap (r) is performed with respect to a cylindrical volume of radius r and height equal to the gap width, d. For clarity of presentation, each curve is shifted vertically by the value indicated in square brackets.
Self-diffusivity as a function of transverse gap position, D(x), for each of the six “wet” confined phases computed at the boxed state points in the phase diagram in Fig. 2. For each phase, the gap region was partitioned into a number of slices of equal width commensurate with the number of lateral layers in the phase and is reflected in the number of points in each curve. The approximate symmetry of the curves around x = 0 nm arises due to the mirror symmetry of the system around the mid-plane of the gap. Uncertainties were computed from the standard deviation in the self-diffusivity computed over the five 4 ns blocks of the full 20 ns trajectory. Where error bars are not visible, the uncertainty is smaller than the symbol size. The inset provides a zoomed in view of the bilayer ice self-diffusivity profile.
Geometry of the monolayer ice. (a) Plan and transverse snapshots of the monolayer ice phase at [P = 0.20 GPa, d = 0.45 nm]. The field of view has been restricted to (−0.25) <x ⩽ 0.25 nm in the x-direction to permit visualization of the inter-plate region. Templating by the walls imparts a hexagonal structure inducing an arrangement of O atoms resembling that of both a single sheet of the bilayer ice and the theoretical structure of a non-dissociated monolayer of water adsorbed onto a Ru(0001) surface.116,117 (b) Since this morphology resides in the plane, satisfaction of all hydrogen bonds is forbidden, prohibiting adherence to the bulk ice rules.110 A simple geometric proof by contradiction of this argument proceeds by assuming that all donor-acceptor pairs are satisfied, which requires that all hydrogen bond around the circumference of a single hexagon point in the same direction, as illustrated in the central cell of this schematic. The inconsistency arises from the imposition upon all adjacent hexagons a configuration of three hydrogen bonds in which two point in one direction and the third in the other. (c) Geometric arguments dictate that at least one defect (unsatisfied hydrogen bond) must exist per hexagonal cell. The number of defects is minimized in an arrangement in which the orientation of hydrogen bond vectors alternates between clockwise (red) and counter-clockwise (blue) along two axes of the hexagonal tiling, generating alternating stripes of a single orientation in the third axis. Defects occur where hexagons of the same color share an edge.
Probability density function for the observation of either H atom of a confined water molecule as a function of the transverse distance across the gap, P H, gap (x), for monolayer ice phase observed at state point [P = 0.20 GPa, d = 0.45 nm]. Normalization by the area of the plates, A = (3.218 × 3.215) nm2 = 10.35 nm2 yields the local density profile across the gap, ρ H, gap (x).
Distribution of O-H vector angles of inclination in the monolayer ice. (a) Probability density function of the O-H vector angle of inclination, θ, with respect to the plane of the gap, averaged over the two O-H vectors of each water molecule residing within the gap, . We adopt the convention in which θ = 90° for an O-H vector pointing towards the upper plate, (–90)° for those pointing towards the lower plate, and 0° for those vectors residing within the y,z-plane. The approximate symmetry around θ = 0° arises from the mirror symmetry of the system around the mid-plane of the gap. (b) Mean field FES experienced by a particular water molecule residing within the monolayer ice, parametrized by the O-H angles of inclination of its O-H1 and O-H2 vectors, , where G is the Gibbs free energy, β = 1/k B T, k B is Boltzmann's constant and T the absolute temperature. The plots in both panels were computed by compiling histograms over the final 4 ns of the equilibrium simulation trajectory at the state point [P = 0.20 GPa, d = 0.45 nm], employing a bin size of Δθ = 2° in (a), and Δθ1 = Δθ2 = 4° in (b). The edges of the landscape in (b) appear blunted due to the cubic grid employed in compiling the histograms.
Nucleation histogram for 117 bilayer ice to vapor decompression-induced sublimation runs that exhibited complete drying transitions. The gap region was discretized into cells of size [Δy = 0.495 nm, Δz = 0.429 nm], resulting in a 6×7 grid in the y,z-plane that spans the gap in the x-direction. Cell size was dictated by the geometry of the cristobalite plates (cf. Fig. 1(d)). The heat map illustrates the number of bridging cavity nucleation events observed within each cell.
The (a) mean binary occupancy vector, , and (b)–(f) top five principal component vectors, , i = 1, …, 5, resulting from application of PCA to the ensemble of 201 decompression-induced sublimation trajectories from the bilayer crystalline solid to vapor at plate separations of d = 0.60 nm. Grid dimensions and orientation are identical to that presented in Fig. 1(d). The scalar value in each cell of these reconstructions may be interpreted as the average probability of that cell being occupied by at least one water molecule O atom, with the understanding that negative values are permissible in panels (b)–(f), since these probabilities do not take on meaning until a linear combination of these PCs and the mean binary occupancy vector is formed. The color bar pertaining to panel (a) is displayed directly to the right of this panel, whereas panels (b)–(f) share the color bar adjacent to panel (f).
Embedding of each snapshot from the 201 sublimation trajectories into the top three PCs of this data set. α i is defined in Eq. (1) as the projection into of the binary occupancy vector associated with a particular snapshot. Each point has been colored according to the fractional occupancy of the gap in that particular snapshot, defined as the fraction of cells in Fig. 1(d) which contain at least one water molecule O atom.
To facilitate visualization of the internal structure of the three-dimensional point cloud in Fig. 12, the point cloud was cut along the axis containing its projection into nine slices over the α1 ranges (a) [–1.718, −0.999), (b) [−0.999, −0.280), (c) [−0.280,0.440] (d) [0.440,1.159), (e) [1.159,1.878), (f) [1.878,2.597), (g) [2.597,3.316), (h) [3.316,4.036), and (i) [4.036,4.755). The two-dimensional probability distributions in this figure were generated by collecting histograms over the PC2-PC3 projection of the points in each slice (Fig. S8107), using cubic mesh with side length 0.5. The color bar indicates the total probability, P, residing in each cell of the mesh. Interpolative contours in P are plotted in increments of 0.01. The × symbol in each frame (a)–(i) marks the location of representative system snapshots presented in the corresponding panels of Fig. 14.
Sequence of nine system snapshots illustrating the progression of drying from the (a) bilayer ice to (i) vapor phase. The plates have been removed from view to better show the gap region. The process commences by the nucleation of a vapor cavity near one corner of the plates, which proceeds to expand down one side, and then another, leaving a promontory of bilayer ice extending into the gap region. Drying completes when this promontory undergoes fluctuation induced melting and recedes into the bulk fluid. In some trajectories, a small island of water molecules transiently exists the center of the gap, as shown in panel (h). These snapshots were selected as representative configurations extracted at the locations indicated by × symbols in the corresponding panels in Fig. 13.
Article metrics loading...
Full text loading...