| 1 |
%load any "/usepackage" here... |
| 2 |
\documentclass{pnastwo} |
| 3 |
|
| 4 |
%% ADDITIONAL OPTIONAL STYLE FILES Font specification |
| 5 |
|
| 6 |
%\usepackage{PNASTWOF} |
| 7 |
\usepackage[version=3]{mhchem} |
| 8 |
\usepackage[round,numbers,sort&compress]{natbib} |
| 9 |
\usepackage{fixltx2e} |
| 10 |
\usepackage{booktabs} |
| 11 |
\usepackage{multirow} |
| 12 |
\usepackage{tablefootnote} |
| 13 |
|
| 14 |
\bibpunct{(}{)}{,}{n}{,}{,} |
| 15 |
\bibliographystyle{pnas2011} |
| 16 |
|
| 17 |
\renewcommand{\thefigure}{S\arabic{figure}} |
| 18 |
\renewcommand{\thetable}{S\arabic{table}} |
| 19 |
\renewcommand{\theequation}{S\arabic{equation}} |
| 20 |
\renewcommand{\thesection}{S\arabic{section}} |
| 21 |
|
| 22 |
%% OPTIONAL MACRO DEFINITIONS |
| 23 |
\def\s{\sigma} |
| 24 |
%%%%%%%%%%%% |
| 25 |
%% For PNAS Only: |
| 26 |
%\url{www.pnas.org/cgi/doi/10.1073/pnas.0709640104} |
| 27 |
\copyrightyear{2014} |
| 28 |
\issuedate{Issue Date} |
| 29 |
\volume{Volume} |
| 30 |
\issuenumber{Issue Number} |
| 31 |
%\setcounter{page}{2687} %Set page number here if desired |
| 32 |
%%%%%%%%%%%% |
| 33 |
|
| 34 |
\begin{document} |
| 35 |
|
| 36 |
\title{Supporting Information for: \\ |
| 37 |
The different facets of ice have different hydrophilicities: Friction at water / |
| 38 |
ice-I\textsubscript{h} interfaces} |
| 39 |
|
| 40 |
\author{Patrick B. Louden\affil{1}{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, |
| 41 |
IN 46556} |
| 42 |
\and |
| 43 |
J. Daniel Gezelter\affil{1}{}} |
| 44 |
|
| 45 |
\contributor{Submitted to Proceedings of the National Academy of Sciences |
| 46 |
of the United States of America} |
| 47 |
|
| 48 |
\maketitle |
| 49 |
|
| 50 |
\begin{article} |
| 51 |
|
| 52 |
\section{Overview} |
| 53 |
The supporting information contains further details about the model |
| 54 |
construction, analysis methods, and supplies figures that support the |
| 55 |
data presented in the main text. |
| 56 |
|
| 57 |
\section{Construction of the Ice / Water interfaces} |
| 58 |
Ice I$_\mathrm{h}$ crystallizes in the hexagonal space group |
| 59 |
P$6_3/mmc$, and common ice crystals form hexagonal plates with the |
| 60 |
basal face, $\{0~0~0~1\}$, forming the top and bottom of each plate, and |
| 61 |
the prismatic facet, $\{1~0~\bar{1}~0\}$, forming the sides. In extreme |
| 62 |
temperatures or low water saturation conditions, ice crystals can |
| 63 |
easily form as hollow columns, needles and dendrites. These are |
| 64 |
structures that expose other crystalline facets of the ice to the |
| 65 |
surroundings. Among the more common facets are the secondary prism, |
| 66 |
$\{1~1~\bar{2}~0\}$, and pyramidal, $\{2~0~\bar{2}~1\}$, faces. |
| 67 |
|
| 68 |
We found it most useful to work with proton-ordered, zero-dipole |
| 69 |
crystals that expose strips of dangling H-atoms and lone |
| 70 |
pairs~\cite{Buch:2008fk}. Our structures were created starting from |
| 71 |
Structure 6 of Hirsch and Ojam\"{a}e's set of orthorhombic |
| 72 |
representations for ice-I$_{h}$~\cite{Hirsch04}. This crystal |
| 73 |
structure was cleaved along the four different facets. The exposed |
| 74 |
face was reoriented normal to the $z$-axis of the simulation cell, and |
| 75 |
the structures were extended to form large exposed facets in |
| 76 |
rectangular box geometries. Liquid water boxes were created with |
| 77 |
identical dimensions (in $x$ and $y$) as the ice, with a $z$ dimension |
| 78 |
of three times that of the ice block, and a density corresponding to 1 |
| 79 |
g / cm$^3$. Each of the ice slabs and water boxes were independently |
| 80 |
equilibrated at a pressure of 1 atm, and the resulting systems were |
| 81 |
merged by carving out any liquid water molecules within 3 \AA\ of any |
| 82 |
atoms in the ice slabs. Each of the combined ice/water systems were |
| 83 |
then equilibrated at 225K, which is the liquid-ice coexistence |
| 84 |
temperature for SPC/E water~\cite{Bryk02}. Reference |
| 85 |
\citealp{Louden13} contains a more detailed explanation of the |
| 86 |
construction of similar ice/water interfaces. The resulting dimensions |
| 87 |
as well as the number of ice and liquid water molecules contained in |
| 88 |
each of these systems are shown in Table S1. |
| 89 |
|
| 90 |
\section{Further details on the shearing (RNEMD) simulations} |
| 91 |
All simulations were performed using OpenMD~\cite{OOPSE,openmd}, with a |
| 92 |
time step of 2 fs and periodic boundary conditions in all three |
| 93 |
dimensions. Electrostatics were handled using the damped-shifted |
| 94 |
force real-space electrostatic kernel~\cite{Ewald}. The systems were |
| 95 |
divided into 100 bins along the $z$-axis for the VSS-RNEMD moves, |
| 96 |
which were attempted every 2~fs. |
| 97 |
|
| 98 |
The interfaces were equilibrated for a total of 10 ns at equilibrium |
| 99 |
conditions before being exposed to either a shear or thermal gradient. |
| 100 |
This consisted of 5 ns under a constant temperature (NVT) integrator |
| 101 |
set to 225~K followed by 5 ns under a microcanonical (NVE) integrator. |
| 102 |
Weak thermal gradients were allowed to develop using the VSS-RNEMD |
| 103 |
(NVE) integrator using a small thermal flux ($-2.0\times 10^{-6}$ |
| 104 |
kcal/mol/\AA$^2$/fs) for a duration of 5 ns to allow the gradient to |
| 105 |
stabilize. The resulting temperature gradient was $\approx$ 10K over |
| 106 |
the entire box length, which was sufficient to keep the temperature at |
| 107 |
the interface within $\pm 1$ K of the 225~K target. |
| 108 |
|
| 109 |
Velocity gradients were then imposed using the VSS-RNEMD (NVE) |
| 110 |
integrator with a range of momentum fluxes. These gradients were |
| 111 |
allowed to stabilize for 1~ns before data collection started. Once |
| 112 |
established, four successive 0.5~ns runs were performed for each shear |
| 113 |
rate. During these simulations, configurations of the system were |
| 114 |
stored every 1~ps, and statistics on the structure and dynamics in |
| 115 |
each bin were accumulated throughout the simulations. Although there |
| 116 |
was some small variation in the measured interfacial width between |
| 117 |
succcessive runs, no indication of bulk melting or crystallization was |
| 118 |
observed. That is, no large scale changes in the positions of the top |
| 119 |
and bottom interfaces occurred during the simulations. |
| 120 |
|
| 121 |
\section{A second method for computing contact angles} |
| 122 |
In addition to the spherical cap method outlined in the main text, a |
| 123 |
second method for obtaining the contact angle was described by |
| 124 |
Ruijter, Blake, and Coninck~\cite{Ruijter99}. This method uses a |
| 125 |
cylindrical averaging of the droplet's density profile. A threshold |
| 126 |
density of 0.5 g cm\textsuperscript{-3} is used to estimate the |
| 127 |
location of the edge of the droplet. The $r$ and $z$-dependence of |
| 128 |
the droplet's edge is then fit to a circle, and the contact angle is |
| 129 |
computed from the intersection of the fit circle with the $z$-axis |
| 130 |
location of the solid surface. Again, for each stored configuration, |
| 131 |
the density profile in a set of annular shells was computed. Due to |
| 132 |
large density fluctuations close to the ice, all shells located within |
| 133 |
2 \AA\ of the ice surface were left out of the circular fits. The |
| 134 |
height of the solid surface ($z_\mathrm{suface}$) along with the best |
| 135 |
fitting origin ($z_\mathrm{droplet}$) and radius |
| 136 |
($r_\mathrm{droplet}$) of the droplet can then be used to compute the |
| 137 |
contact angle, |
| 138 |
\begin{equation} |
| 139 |
\theta = 90 + \frac{180}{\pi} \sin^{-1}\left(\frac{z_\mathrm{droplet} - |
| 140 |
z_\mathrm{surface}}{r_\mathrm{droplet}} \right). |
| 141 |
\end{equation} |
| 142 |
|
| 143 |
\section{Interfacial widths using structural information} |
| 144 |
To determine the structural widths of the interfaces under shear, each |
| 145 |
of the systems was divided into 100 bins along the $z$-dimension, and |
| 146 |
the local tetrahedral order parameter (Eq. 5 in Reference |
| 147 |
\citealp{Louden13}) was time-averaged in each bin for the duration of |
| 148 |
the shearing simulation. The spatial dependence of this order |
| 149 |
parameter, $q(z)$, is the tetrahedrality profile of the interface. |
| 150 |
The lower panels in figures S2-S5 in the SI show tetrahedrality |
| 151 |
profiles (in circles) for each of the four interfaces. The $q(z)$ |
| 152 |
function has a range of $(0,1)$, where a value of unity indicates a |
| 153 |
perfectly tetrahedral environment. The $q(z)$ for the bulk liquid was |
| 154 |
found to be $\approx~0.77$, while values of $\approx~0.92$ were more |
| 155 |
common in the ice. The tetrahedrality profiles were fit using a |
| 156 |
hyperbolic tangent function (see Eq. 6 in Reference |
| 157 |
\citealp{Louden13}) designed to smoothly fit the bulk to ice |
| 158 |
transition while accounting for the weak thermal gradient. In panels |
| 159 |
$b$ and $c$ of the same figures, the resulting thermal and velocity |
| 160 |
gradients from an imposed kinetic energy and momentum fluxes can be |
| 161 |
seen. The vertical dotted lines traversing these figures indicate the |
| 162 |
midpoints of the interfaces as determined by the tetrahedrality |
| 163 |
profiles. |
| 164 |
|
| 165 |
\section{Interfacial widths using dynamic information} |
| 166 |
To determine the dynamic widths of the interfaces under shear, each of |
| 167 |
the systems was divided into bins along the $z$-dimension ($\approx$ 3 |
| 168 |
\AA\ wide) and $C_2(z,t)$ was computed using only those molecules that |
| 169 |
were in the bin at the initial time. To compute these correlation |
| 170 |
functions, each of the 0.5 ns simulations was followed by a shorter |
| 171 |
200 ps microcanonical (NVE) simulation in which the positions and |
| 172 |
orientations of every molecule in the system were recorded every 0.1 |
| 173 |
ps. |
| 174 |
|
| 175 |
The time-dependence was fit to a triexponential decay, with three time |
| 176 |
constants: $\tau_{short}$, measuring the librational motion of the |
| 177 |
water molecules, $\tau_{middle}$, measuring the timescale for breaking |
| 178 |
and making of hydrogen bonds, and $\tau_{long}$, corresponding to the |
| 179 |
translational motion of the water molecules. An additional constant |
| 180 |
was introduced in the fits to describe molecules in the crystal which |
| 181 |
do not experience long-time orientational decay. |
| 182 |
|
| 183 |
In Figures S6-S9, the $z$-coordinate profiles for the three decay |
| 184 |
constants, $\tau_{short}$, $\tau_{middle}$, and $\tau_{long}$ for the |
| 185 |
different interfaces are shown. (Figures S6 \& S7 are new results, |
| 186 |
and Figures S8 \& S9 are updated plots from Ref \citealp{Louden13}.) |
| 187 |
In the liquid regions of all four interfaces, we observe |
| 188 |
$\tau_{middle}$ and $\tau_{long}$ to have approximately consistent |
| 189 |
values of $3-6$ ps and $30-40$ ps, respectively. Both of these times |
| 190 |
increase in value approaching the interface. Approaching the |
| 191 |
interface, we also observe that $\tau_{short}$ decreases from its |
| 192 |
liquid-state value of $72-76$ fs. The approximate values for the |
| 193 |
decay constants and the trends approaching the interface match those |
| 194 |
reported previously for the basal and prismatic interfaces. |
| 195 |
|
| 196 |
We have estimated the dynamic interfacial width $d_\mathrm{dyn}$ by |
| 197 |
fitting the profiles of all the three orientational time constants |
| 198 |
with an exponential decay to the bulk-liquid behavior, |
| 199 |
\begin{equation}\label{tauFit} |
| 200 |
\tau(z)\approx\tau_{liquid}+(\tau_{wall}-\tau_{liquid})e^{-(z-z_{wall})/d_\mathrm{dyn}} |
| 201 |
\end{equation} |
| 202 |
where $\tau_{liquid}$ and $\tau_{wall}$ are the liquid and projected |
| 203 |
wall values of the decay constants, $z_{wall}$ is the location of the |
| 204 |
interface, as measured by the structural order parameter. These |
| 205 |
values are shown in table 1 in the main text. Because the bins must be |
| 206 |
quite wide to obtain reasonable profiles of $C_2(z,t)$, the error |
| 207 |
estimates for the dynamic widths of the interface are significantly |
| 208 |
larger than for the structural widths. However, all four interfaces |
| 209 |
exhibit dynamic widths that are significantly below 1~nm, and are in |
| 210 |
reasonable agreement with the structural width above. |
| 211 |
|
| 212 |
\bibliography{iceWater} |
| 213 |
\end{article} |
| 214 |
|
| 215 |
\begin{table}[h] |
| 216 |
\centering |
| 217 |
\caption{Sizes of the droplet and shearing simulations. Cell |
| 218 |
dimensions are measured in \AA. \label{tab:method}} |
| 219 |
\begin{tabular}{r|cccc|ccccc} |
| 220 |
\toprule |
| 221 |
\multirow{2}{*}{Interface} & \multicolumn{4}{c|}{Droplet} & \multicolumn{5}{c}{Shearing} \\ |
| 222 |
& $N_\mathrm{ice}$ & $N_\mathrm{droplet}$ & $L_x$ & $L_y$ & $N_\mathrm{ice}$ & $N_\mathrm{liquid}$ & $L_x$ & $L_y$ & $L_z$ \\ |
| 223 |
\midrule |
| 224 |
Basal $\{0001\}$ & 12960 & 2048 & 134.70 & 140.04 & 900 & 1846 & 23.87 & 35.83 & 98.64 \\ |
| 225 |
Pyramidal $\{2~0~\bar{2}~1\}$ & 11136 & 2048 & 143.75 & 121.41 & 1216 & 2203 & 37.47 & 29.50 & 93.02 \\ |
| 226 |
Prismatic $\{1~0~\bar{1}~0\}$ & 9900 & 2048 & 110.04 & 115.00 & 3000 & 5464 & 35.95 & 35.65 & 205.77 \\ |
| 227 |
Secondary Prism $\{1~1~\bar{2}~0\}$ & 11520 & 2048 & 146.72 & 124.48 & 3840 & 8176 & 71.87 & 31.66 & 161.55 \\ |
| 228 |
\bottomrule |
| 229 |
\end{tabular} |
| 230 |
\end{table} |
| 231 |
|
| 232 |
%S1: contact angle |
| 233 |
\begin{figure} |
| 234 |
\includegraphics[width=\linewidth]{ContactAngle} |
| 235 |
\caption{\label{fig:ContactAngle} The dynamic contact angle of a |
| 236 |
droplet after approaching each of the four ice facets. The decay to |
| 237 |
an equilibrium contact angle displays similar dynamics. Although |
| 238 |
all the surfaces are hydrophilic, the long-time behavior stabilizes |
| 239 |
to significantly flatter droplets for the basal and pyramidal |
| 240 |
facets. This suggests a difference in hydrophilicity for these |
| 241 |
facets compared with the two prismatic facets.} |
| 242 |
\end{figure} |
| 243 |
|
| 244 |
|
| 245 |
%S2-S5 are the z-rnemd profiles |
| 246 |
\begin{figure} |
| 247 |
\includegraphics[width=\linewidth]{Pyr_comic_strip} |
| 248 |
\caption{\label{fig:pyrComic} Properties of the pyramidal interface |
| 249 |
being sheared through water at 3.8 ms\textsuperscript{-1}. Lower |
| 250 |
panel: the local tetrahedral order parameter, $q(z)$, (circles) and |
| 251 |
the hyperbolic tangent fit (turquoise line). Middle panel: the |
| 252 |
imposed thermal gradient required to maintain a fixed interfacial |
| 253 |
temperature of 225 K. Upper panel: the transverse velocity gradient |
| 254 |
that develops in response to an imposed momentum flux. The vertical |
| 255 |
dotted lines indicate the locations of the midpoints of the two |
| 256 |
interfaces.} |
| 257 |
\end{figure} |
| 258 |
|
| 259 |
\begin{figure} |
| 260 |
\includegraphics[width=\linewidth]{SP_comic_strip} |
| 261 |
\caption{\label{fig:spComic} The secondary prism interface with a shear |
| 262 |
rate of 3.5 \ |
| 263 |
ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:pyrComic}.} |
| 264 |
\end{figure} |
| 265 |
|
| 266 |
\begin{figure} |
| 267 |
\includegraphics[width=\linewidth]{B_comic_strip} |
| 268 |
\caption{\label{fig:bComic} The basal interface with a shear |
| 269 |
rate of 1.3 \ |
| 270 |
ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:pyrComic}.} |
| 271 |
\end{figure} |
| 272 |
|
| 273 |
\begin{figure} |
| 274 |
\includegraphics[width=\linewidth]{prismatic_comic_strip} |
| 275 |
\caption{\label{fig:pComic} The prismatic interface with a shear |
| 276 |
rate of 2 \ |
| 277 |
ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:pyrComic}.} |
| 278 |
\end{figure} |
| 279 |
|
| 280 |
%Figures S6-S9 are the z-orientation times |
| 281 |
\begin{figure} |
| 282 |
\includegraphics[width=\linewidth]{Pyr-orient} |
| 283 |
\caption{\label{fig:PyrOrient} The three decay constants of the |
| 284 |
orientational time correlation function, $C_2(z,t)$, for water as a |
| 285 |
function of distance from the center of the ice slab. The vertical |
| 286 |
dashed line indicates the edge of the pyramidal ice slab determined |
| 287 |
by the local order tetrahedral parameter. The control (circles) and |
| 288 |
sheared (squares) simulations were fit using shifted-exponential |
| 289 |
decay (see Eq. 9 in Ref. \citealp{Louden13}).} |
| 290 |
\end{figure} |
| 291 |
|
| 292 |
\begin{figure} |
| 293 |
\includegraphics[width=\linewidth]{SP-orient} |
| 294 |
\caption{\label{fig:SPorient} Decay constants for $C_2(z,t)$ at the secondary |
| 295 |
prism face. Panel descriptions match those in \ref{fig:PyrOrient}.} |
| 296 |
\end{figure} |
| 297 |
|
| 298 |
|
| 299 |
\begin{figure} |
| 300 |
\includegraphics[width=\linewidth]{B-orient} |
| 301 |
\caption{\label{fig:Borient} Decay constants for $C_2(z,t)$ at the basal face. Panel descriptions match those in \ref{fig:PyrOrient}.} |
| 302 |
\end{figure} |
| 303 |
|
| 304 |
\begin{figure} |
| 305 |
\includegraphics[width=\linewidth]{prismatic-orient} |
| 306 |
\caption{\label{fig:Porient} Decay constants for $C_2(z,t)$ at the |
| 307 |
prismatic face. Panel descriptions match those in \ref{fig:PyrOrient}.} |
| 308 |
\end{figure} |
| 309 |
|
| 310 |
|
| 311 |
\end{document} |
