1 |
gezelter |
3914 |
\documentclass[journal = jpccck, manuscript = article]{achemso} |
2 |
|
|
\setkeys{acs}{usetitle = true} |
3 |
|
|
\usepackage{achemso} |
4 |
|
|
\usepackage{natbib} |
5 |
|
|
\usepackage{multirow} |
6 |
|
|
\usepackage{wrapfig} |
7 |
|
|
\usepackage{fixltx2e} |
8 |
|
|
%\mciteErrorOnUnknownfalse |
9 |
gezelter |
3897 |
|
10 |
gezelter |
3914 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
11 |
|
|
\usepackage{url} |
12 |
gezelter |
3897 |
|
13 |
gezelter |
3949 |
\title{Simulations of solid-liquid friction at ice-I$_\mathrm{h}$ / |
14 |
|
|
water interfaces} |
15 |
gezelter |
3897 |
|
16 |
|
|
\author{P. B. Louden} |
17 |
gezelter |
3914 |
\author{J. Daniel Gezelter} |
18 |
|
|
\email{gezelter@nd.edu} |
19 |
|
|
\affiliation[University of Notre Dame]{251 Nieuwland Science Hall\\ |
20 |
|
|
Department of Chemistry and Biochemistry\\ University of Notre |
21 |
|
|
Dame\\ Notre Dame, Indiana 46556} |
22 |
gezelter |
3897 |
|
23 |
gezelter |
3914 |
\keywords{} |
24 |
gezelter |
3897 |
|
25 |
gezelter |
3914 |
\begin{document} |
26 |
gezelter |
3897 |
|
27 |
gezelter |
3914 |
\begin{abstract} |
28 |
gezelter |
3945 |
We have investigated the structural and dynamic properties of the |
29 |
gezelter |
3954 |
basal and prismatic facets of the ice I$_\mathrm{h}$ / water |
30 |
|
|
interface when the solid phase is drawn through the liquid |
31 |
gezelter |
3945 |
(i.e. sheared relative to the fluid phase). To impose the shear, we |
32 |
|
|
utilized a velocity-shearing and scaling (VSS) approach to reverse |
33 |
|
|
non-equilibrium molecular dynamics (RNEMD). This method can create |
34 |
|
|
simultaneous temperature and velocity gradients and allow the |
35 |
|
|
measurement of transport properties at interfaces. The interfacial |
36 |
gezelter |
3954 |
width was found to be independent of the relative velocity of the |
37 |
|
|
ice and liquid layers over a wide range of shear rates. Decays of |
38 |
|
|
molecular orientational time correlation functions gave similar |
39 |
|
|
estimates for the width of the interfaces, although the short- and |
40 |
|
|
longer-time decay components behave differently closer to the |
41 |
|
|
interface. Although both facets of ice are in ``stick'' boundary |
42 |
|
|
conditions in liquid water, the solid-liquid friction coefficients |
43 |
|
|
were found to be significantly different for the basal and prismatic |
44 |
|
|
facets of ice. |
45 |
gezelter |
3914 |
\end{abstract} |
46 |
gezelter |
3897 |
|
47 |
gezelter |
3914 |
\newpage |
48 |
gezelter |
3897 |
|
49 |
|
|
\section{Introduction} |
50 |
plouden |
3919 |
%-----Outline of Intro--------------- |
51 |
|
|
% in general, ice/water interface is important b/c .... |
52 |
|
|
% here are some people who have worked on ice/water, trying to understand the processes above .... |
53 |
|
|
% with the recent development of VSS-RNEMD, we can now look at the shearing problem |
54 |
|
|
% talk about what we will present in this paper |
55 |
|
|
% -------End Intro------------------ |
56 |
gezelter |
3897 |
|
57 |
plouden |
3919 |
%Gay02: cites many other ice/water papers, make sure to cite them. |
58 |
|
|
|
59 |
gezelter |
3945 |
Understanding the ice/water interface is essential for explaining |
60 |
|
|
complex processes such as nucleation and crystal |
61 |
|
|
growth,\cite{Han92,Granasy95,Vanfleet95} crystal |
62 |
|
|
melting,\cite{Weber83,Han92,Sakai96,Sakai96B} and some fascinating |
63 |
|
|
biological processes, such as the behavior of the antifreeze proteins |
64 |
|
|
found in winter flounder,\cite{Wierzbicki07, Chapsky97} and certain |
65 |
|
|
terrestrial arthropods.\cite{Duman:2001qy,Meister29012013} There has |
66 |
|
|
been significant progress on understanding the structure and dynamics |
67 |
|
|
of quiescent ice/water interfaces utilizing both theory and |
68 |
|
|
experiment. Haymet \emph{et al.} have done extensive work on ice I$_\mathrm{h}$, |
69 |
|
|
including characterizing and determining the width of the ice/water |
70 |
plouden |
3950 |
interface for the SPC\cite{Karim90}, SPC/E\cite{Gay02,Bryk02}, CF1\cite{Hayward01,Hayward02}, and TIP4P\cite{Karim88} models for |
71 |
|
|
water. |
72 |
gezelter |
3945 |
More recently, Haymet \emph{et al.} have investigated the effects |
73 |
|
|
cations and anions have on crystal |
74 |
|
|
nucleation\cite{Bryk04,Smith05,Wilson08,Wilson10}. Nada \emph{et al.} |
75 |
|
|
have also studied ice/water |
76 |
|
|
interfaces,\cite{Nada95,Nada00,Nada03,Nada12} and have found that the |
77 |
|
|
differential growth rates of the facets of ice I$_\mathrm{h}$ are due to the the |
78 |
|
|
reordering of the hydrogen bonding network\cite{Nada05}. |
79 |
plouden |
3919 |
|
80 |
gezelter |
3945 |
The movement of liquid water over the facets of ice has been less |
81 |
|
|
thoroughly studied than the quiescent surfaces. This process is |
82 |
|
|
potentially important in understanding transport of large blocks of |
83 |
|
|
ice in water (which has important implications in the earth sciences), |
84 |
|
|
as well as the relative motion of crystal-crystal interfaces that have |
85 |
|
|
been separated by nanometer-scale fluid domains. In addition to |
86 |
|
|
understanding both the structure and thickness of the interfacial |
87 |
|
|
regions, it is important to understand the molecular origin of |
88 |
|
|
friction, drag, and other changes in dynamical properties of the |
89 |
|
|
liquid in the regions close to the surface that are altered by the |
90 |
|
|
presence of a shearing of the bulk fluid relative to the solid phase. |
91 |
plouden |
3919 |
|
92 |
gezelter |
3945 |
In this work, we apply a recently-developed velocity shearing and |
93 |
|
|
scaling approach to reverse non-equilibrium molecular dynamics |
94 |
|
|
(VSS-RNEMD). This method makes it possible to calculate transport |
95 |
|
|
properties like the interfacial thermal conductance across |
96 |
|
|
heterogeneous interfaces,\cite{Kuang12} and can create simultaneous |
97 |
|
|
temperature and velocity gradients and allow the measurement of |
98 |
|
|
friction and thermal transport properties at interfaces. This has |
99 |
|
|
allowed us to investigate the width of the ice/water interface as the |
100 |
|
|
ice is sheared through the liquid, while simultaneously imposing a |
101 |
|
|
weak thermal gradient to prevent frictional heating of the interface. |
102 |
|
|
In the sections that follow, we discuss the methodology for creating |
103 |
|
|
and simulating ice/water interfaces under shear and provide results |
104 |
|
|
from both structural and dynamical correlation functions. We also |
105 |
|
|
show that the solid-liquid interfacial friction coefficient depends |
106 |
|
|
sensitively on the details of the surface morphology. |
107 |
|
|
|
108 |
gezelter |
3897 |
\section{Methodology} |
109 |
|
|
|
110 |
gezelter |
3945 |
\subsection{Stable ice I$_\mathrm{h}$ / water interfaces under shear} |
111 |
gezelter |
3914 |
|
112 |
gezelter |
3946 |
The structure of ice I$_\mathrm{h}$ is well understood; it |
113 |
gezelter |
3914 |
crystallizes in a hexagonal space group P$6_3/mmc$, and the hexagonal |
114 |
|
|
crystals of ice have two faces that are commonly exposed, the basal |
115 |
|
|
face $\{0~0~0~1\}$, which forms the top and bottom of each hexagonal |
116 |
|
|
plate, and the prismatic $\{1~0~\bar{1}~0\}$ face which forms the |
117 |
|
|
sides of the plate. Other less-common, but still important, faces of |
118 |
gezelter |
3945 |
ice I$_\mathrm{h}$ are the secondary prism, $\{1~1~\bar{2}~0\}$, and |
119 |
gezelter |
3946 |
pyramidal, $\{2~0~\bar{2}~1\}$, faces. Ice I$_\mathrm{h}$ is normally |
120 |
|
|
proton disordered in bulk crystals, although the surfaces probably |
121 |
|
|
have a preference for proton ordering along strips of dangling H-atoms |
122 |
|
|
and Oxygen lone pairs.\cite{Buch:2008fk} |
123 |
gezelter |
3914 |
|
124 |
|
|
\begin{wraptable}{r}{3.5in} |
125 |
gezelter |
3945 |
\caption{Mapping between the Miller indices of four facets of ice in |
126 |
|
|
the $P6_3/mmc$ crystal system to the orthorhombic $P2_12_12_1$ |
127 |
|
|
system in reference \protect\cite{Hirsch04}} |
128 |
|
|
\label{tab:equiv} |
129 |
gezelter |
3914 |
\begin{tabular}{|ccc|} \hline |
130 |
|
|
& hexagonal & orthorhombic \\ |
131 |
|
|
& ($P6_3/mmc$) & ($P2_12_12_1$) \\ |
132 |
|
|
crystal face & Miller indices & equivalent \\ \hline |
133 |
|
|
basal & $\{0~0~0~1\}$ & $\{0~0~1\}$ \\ |
134 |
|
|
prism & $\{1~0~\bar{1}~0\}$ & $\{1~0~0\}$ \\ |
135 |
|
|
secondary prism & $\{1~1~\bar{2}~0\}$ & $\{1~3~0\}$ \\ |
136 |
gezelter |
3946 |
pyramidal & $\{2~0~\bar{2}~1\}$ & $\{2~0~1\}$ \\ \hline |
137 |
gezelter |
3914 |
\end{tabular} |
138 |
|
|
\end{wraptable} |
139 |
|
|
|
140 |
gezelter |
3946 |
For small simulated ice interfaces, it is useful to work with |
141 |
|
|
proton-ordered, but zero-dipole crystal that exposes these strips of |
142 |
|
|
dangling H-atoms and lone pairs. When placing another material in |
143 |
|
|
contact with one of the ice crystalline planes, it is also quite |
144 |
|
|
useful to have an orthorhombic (rectangular) box. Recent work by |
145 |
|
|
Hirsch and Ojam\"{a}e describes a number of alternative crystal |
146 |
|
|
systems for proton-ordered bulk ice I$_\mathrm{h}$ using orthorhombic |
147 |
|
|
cells.\cite{Hirsch04} |
148 |
|
|
|
149 |
|
|
In this work, we are using Hirsch and Ojam\"{a}e's structure 6 which |
150 |
|
|
is an orthorhombic cell ($P2_12_12_1$) that produces a proton-ordered |
151 |
|
|
version of ice I$_\mathrm{h}$. Table \ref{tab:equiv} contains a |
152 |
|
|
mapping between the Miller indices of common ice facets in the |
153 |
|
|
P$6_3/mmc$ crystal system and those in the Hirsch and Ojam\"{a}e |
154 |
|
|
$P2_12_12_1$ system. |
155 |
|
|
|
156 |
gezelter |
3914 |
Structure 6 from the Hirsch and Ojam\"{a}e paper has lattice |
157 |
|
|
parameters $a = 4.49225$, $b = 7.78080$, $c = 7.33581$ and two water |
158 |
gezelter |
3945 |
molecules whose atoms reside at fractional coordinates given in table |
159 |
|
|
\ref{tab:p212121}. To construct the basal and prismatic interfaces, |
160 |
|
|
these crystallographic coordinates were used to construct an |
161 |
|
|
orthorhombic unit cell which was then replicated in all three |
162 |
|
|
dimensions yielding a proton-ordered block of ice I$_{h}$. To expose |
163 |
|
|
the desired face, the orthorhombic representation was then cut along |
164 |
|
|
the ($001$) or ($100$) planes for the basal and prismatic faces |
165 |
|
|
respectively. The resulting block was rotated so that the exposed |
166 |
gezelter |
3946 |
faces were aligned with the $z$-axis normal to the exposed face. The |
167 |
|
|
block was then cut along two perpendicular directions in a way that |
168 |
|
|
allowed for perfect periodic replication in the $x$ and $y$ axes, |
169 |
gezelter |
3945 |
creating a slab with either the basal or prismatic faces exposed along |
170 |
|
|
the $z$ axis. The slab was then replicated in the $x$ and $y$ |
171 |
|
|
dimensions until a desired sample size was obtained. |
172 |
gezelter |
3914 |
|
173 |
gezelter |
3946 |
\begin{wraptable}{r}{2.85in} |
174 |
gezelter |
3945 |
\caption{Fractional coordinates for water in the orthorhombic |
175 |
|
|
$P2_12_12_1$ system for ice I$_\mathrm{h}$ in reference |
176 |
|
|
\protect\cite{Hirsch04}} |
177 |
|
|
\label{tab:p212121} |
178 |
|
|
\begin{tabular}{|cccc|} \hline |
179 |
|
|
atom type & x & y & z \\ \hline |
180 |
|
|
O & 0.75 & 0.1667 & 0.4375 \\ |
181 |
|
|
H & 0.5735 & 0.2202 & 0.4836 \\ |
182 |
|
|
H & 0.7420 & 0.0517 & 0.4836 \\ |
183 |
|
|
O & 0.25 & 0.6667 & 0.4375 \\ |
184 |
|
|
H & 0.2580 & 0.6693 & 0.3071 \\ |
185 |
|
|
H & 0.4265 & 0.7255 & 0.4756 \\ \hline |
186 |
gezelter |
3914 |
\end{tabular} |
187 |
|
|
\end{wraptable} |
188 |
|
|
|
189 |
gezelter |
3946 |
Our ice / water interfaces were created using a box of liquid water |
190 |
|
|
that had the same dimensions (in $x$ and $y$) as the ice block. |
191 |
|
|
Although the experimental solid/liquid coexistence temperature under |
192 |
|
|
atmospheric pressure is close to 273K, Haymet \emph{et al.} have done |
193 |
|
|
extensive work on characterizing the ice/water interface, and find |
194 |
|
|
that the coexistence temperature for simulated water is often quite a |
195 |
|
|
bit different.\cite{Karim88,Karim90,Hayward01,Bryk02,Hayward02} They |
196 |
|
|
have found that for the SPC/E water model,\cite{Berendsen87} which is |
197 |
|
|
also used in this study, the ice/water interface is most stable at |
198 |
|
|
225$\pm$5K.\cite{Bryk02} This liquid box was therefore equilibrated at |
199 |
|
|
225 K and 1 atm of pressure in the NPAT ensemble (with the $z$ axis |
200 |
|
|
allowed to fluctuate to equilibrate to the correct pressure). The |
201 |
|
|
liquid and solid systems were combined by carving out any water |
202 |
|
|
molecule from the liquid simulation cell that was within 3 \AA\ of any |
203 |
|
|
atom in the ice slab. |
204 |
gezelter |
3914 |
|
205 |
|
|
Molecular translation and orientational restraints were applied in the |
206 |
|
|
early stages of equilibration to prevent melting of the ice slab. |
207 |
|
|
These restraints were removed during NVT equilibration, well before |
208 |
gezelter |
3945 |
data collection was carried out. |
209 |
gezelter |
3914 |
|
210 |
gezelter |
3918 |
\subsection{Shearing ice / water interfaces without bulk melting} |
211 |
|
|
|
212 |
gezelter |
3945 |
As a solid is dragged through a liquid, there is frictional heating |
213 |
|
|
that will act to melt the interface. To study the behavior of the |
214 |
|
|
interface under a shear stress without causing the interface to melt, |
215 |
gezelter |
3946 |
it is necessary to apply a weak thermal gradient in combination with |
216 |
|
|
the momentum gradient. This can be accomplished using the velocity |
217 |
gezelter |
3945 |
shearing and scaling (VSS) variant of reverse non-equilibrium |
218 |
|
|
molecular dynamics (RNEMD), which utilizes a series of simultaneous |
219 |
|
|
velocity exchanges between two regions within the simulation |
220 |
|
|
cell.\cite{Kuang12} One of these regions is centered within the ice |
221 |
gezelter |
3946 |
slab, while the other is centrally located in the liquid |
222 |
gezelter |
3918 |
region. VSS-RNEMD provides a set of conservation constraints for |
223 |
gezelter |
3946 |
creating either a momentum flux or a thermal flux (or both |
224 |
|
|
simultaneously) between the two slabs. Satisfying the constraint |
225 |
|
|
equations ensures that the new configurations are sampled from the |
226 |
|
|
same NVE ensemble as before the VSS move. |
227 |
gezelter |
3918 |
|
228 |
|
|
The VSS moves are applied periodically to scale and shift the particle |
229 |
|
|
velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and |
230 |
|
|
$C$) which are separated by half of the simulation box, |
231 |
|
|
\begin{displaymath} |
232 |
|
|
\begin{array}{rclcl} |
233 |
|
|
|
234 |
|
|
& \underline{\mathrm{shearing}} & & |
235 |
|
|
\underline{~~~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~~~} \\ |
236 |
|
|
\mathbf{v}_i \leftarrow & |
237 |
|
|
\mathbf{a}_c\ & + & c\cdot\left(\mathbf{v}_i - \langle\mathbf{v}_c |
238 |
|
|
\rangle\right) + \langle\mathbf{v}_c\rangle \\ |
239 |
|
|
\mathbf{v}_j \leftarrow & |
240 |
|
|
\mathbf{a}_h & + & h\cdot\left(\mathbf{v}_j - \langle\mathbf{v}_h |
241 |
|
|
\rangle\right) + \langle\mathbf{v}_h\rangle . |
242 |
|
|
|
243 |
|
|
\end{array} |
244 |
|
|
\end{displaymath} |
245 |
|
|
Here $\langle\mathbf{v}_c\rangle$ and $\langle\mathbf{v}_h\rangle$ are |
246 |
|
|
the center of mass velocities in the $C$ and $H$ slabs, respectively. |
247 |
|
|
Within the two slabs, particles receive incremental changes or a |
248 |
|
|
``shear'' to their velocities. The amount of shear is governed by the |
249 |
|
|
imposed momentum flux, $\mathbf{j}_z(\mathbf{p})$ |
250 |
|
|
\begin{eqnarray} |
251 |
|
|
\mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\ |
252 |
|
|
\mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2} |
253 |
|
|
\end{eqnarray} |
254 |
|
|
where $M_{\{c,h\}}$ is the total mass of particles within each of the |
255 |
|
|
slabs and $\Delta t$ is the interval between two separate operations. |
256 |
|
|
|
257 |
|
|
To simultaneously impose a thermal flux ($J_z$) between the slabs we |
258 |
|
|
use energy conservation constraints, |
259 |
|
|
\begin{eqnarray} |
260 |
|
|
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\mathbf{v}_c |
261 |
|
|
\rangle^2) + \frac{1}{2}M_c (\langle \mathbf{v}_c \rangle + \mathbf{a}_c)^2 \label{vss3}\\ |
262 |
|
|
K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\mathbf{v}_h |
263 |
|
|
\rangle^2) + \frac{1}{2}M_h (\langle \mathbf{v}_h \rangle + |
264 |
|
|
\mathbf{a}_h)^2 \label{vss4}. |
265 |
|
|
\label{constraint} |
266 |
|
|
\end{eqnarray} |
267 |
|
|
Simultaneous solution of these quadratic formulae for the scaling |
268 |
|
|
coefficients, $c$ and $h$, will ensure that the simulation samples from |
269 |
|
|
the original microcanonical (NVE) ensemble. Here $K_{\{c,h\}}$ is the |
270 |
|
|
instantaneous translational kinetic energy of each slab. At each time |
271 |
|
|
interval, it is a simple matter to solve for $c$, $h$, $\mathbf{a}_c$, |
272 |
|
|
and $\mathbf{a}_h$, subject to the imposed momentum flux, |
273 |
|
|
$j_z(\mathbf{p})$, and thermal flux, $J_z$, values. Since the VSS |
274 |
|
|
operations do not change the kinetic energy due to orientational |
275 |
|
|
degrees of freedom or the potential energy of a system, configurations |
276 |
|
|
after the VSS move have exactly the same energy (and linear |
277 |
|
|
momentum) as before the move. |
278 |
|
|
|
279 |
|
|
As the simulation progresses, the VSS moves are performed on a regular |
280 |
|
|
basis, and the system develops a thermal and/or velocity gradient in |
281 |
gezelter |
3945 |
response to the applied flux. In a bulk material, it is quite simple |
282 |
gezelter |
3918 |
to use the slope of the temperature or velocity gradients to obtain |
283 |
gezelter |
3945 |
either the thermal conductivity or shear viscosity. |
284 |
gezelter |
3918 |
|
285 |
|
|
The VSS-RNEMD approach is versatile in that it may be used to |
286 |
|
|
implement thermal and shear transport simultaneously. Perturbations |
287 |
|
|
of velocities away from the ideal Maxwell-Boltzmann distributions are |
288 |
|
|
minimal, as is thermal anisotropy. This ability to generate |
289 |
|
|
simultaneous thermal and shear fluxes has been previously utilized to |
290 |
|
|
map out the shear viscosity of SPC/E water over a wide range of |
291 |
|
|
temperatures (90~K) with a single 1 ns simulation.\cite{Kuang12} |
292 |
|
|
|
293 |
gezelter |
3945 |
For this work, we are using the VSS-RNEMD method primarily to generate |
294 |
|
|
a shear between the ice slab and the liquid phase, while using a weak |
295 |
gezelter |
3946 |
thermal gradient to maintain the interface at the 225K target |
296 |
|
|
value. This ensures minimal melting of the bulk ice phase and allows |
297 |
|
|
us to control the exact temperature of the interface. |
298 |
gezelter |
3918 |
|
299 |
gezelter |
3897 |
\subsection{Computational Details} |
300 |
gezelter |
3945 |
All simulations were performed using OpenMD,\cite{OOPSE,openmd} with a |
301 |
|
|
time step of 2 fs and periodic boundary conditions in all three |
302 |
|
|
dimensions. Electrostatics were handled using the damped-shifted |
303 |
|
|
force real-space electrostatic kernel.\cite{Ewald} The systems were |
304 |
|
|
divided into 100 bins along the $z$-axis for the VSS-RNEMD moves, |
305 |
|
|
which were attempted every 50 fs. |
306 |
gezelter |
3897 |
|
307 |
gezelter |
3945 |
The interfaces were equilibrated for a total of 10 ns at equilibrium |
308 |
|
|
conditions before being exposed to either a shear or thermal gradient. |
309 |
|
|
This consisted of 5 ns under a constant temperature (NVT) integrator |
310 |
|
|
set to 225K followed by 5 ns under a microcanonical integrator. Weak |
311 |
|
|
thermal gradients were allowed to develop using the VSS-RNEMD (NVE) |
312 |
|
|
integrator using a a small thermal flux ($-2.0\times 10^{-6}$ |
313 |
|
|
kcal/mol/\AA$^2$/fs) for a duration of 5 ns to allow the gradient to |
314 |
gezelter |
3946 |
stabilize. The resulting temperature gradient was $\approx$ 10K over |
315 |
|
|
the entire 100 \AA\ box length, which was sufficient to keep the |
316 |
gezelter |
3945 |
temperature at the interface within $\pm 1$ K of the 225K target. |
317 |
gezelter |
3897 |
|
318 |
gezelter |
3945 |
Velocity gradients were then imposed using the VSS-RNEMD (NVE) |
319 |
|
|
integrator with a range of momentum fluxes. These gradients were |
320 |
|
|
allowed to stabilize for 1 ns before data collection began. Once |
321 |
|
|
established, four successive 0.5 ns runs were performed for each shear |
322 |
|
|
rate. During these simulations, snapshots of the system were taken |
323 |
|
|
every 1 ps, and statistics on the structure and dynamics in each bin |
324 |
|
|
were accumulated throughout the simulations. |
325 |
|
|
|
326 |
plouden |
3941 |
\section{Results and discussion} |
327 |
|
|
|
328 |
gezelter |
3946 |
\subsection{Interfacial width} |
329 |
|
|
Any order parameter or time correlation function that changes as one |
330 |
|
|
crosses an interface from a bulk liquid to a solid can be used to |
331 |
|
|
measure the width of the interface. In previous work on the ice/water |
332 |
gezelter |
3954 |
interface, Haymet {\it et al.}\cite{Bryk02} have utilized structural |
333 |
gezelter |
3946 |
features (including the density) as well as dynamic properties |
334 |
|
|
(including the diffusion constant) to estimate the width of the |
335 |
|
|
interfaces for a number of facets of the ice crystals. Because |
336 |
|
|
VSS-RNEMD imposes a lateral flow, parameters that depend on |
337 |
|
|
translational motion of the molecules (e.g. diffusion) may be |
338 |
gezelter |
3954 |
artificially skewed by the RNEMD moves. A structural parameter is not |
339 |
gezelter |
3946 |
influenced by the RNEMD perturbations to the same degree. Here, we |
340 |
gezelter |
3954 |
have used the local tetrahedral order parameter as described by |
341 |
gezelter |
3946 |
Kumar\cite{Kumar09} and Errington\cite{Errington01} as our principal |
342 |
gezelter |
3945 |
measure of the interfacial width. |
343 |
plouden |
3936 |
|
344 |
gezelter |
3945 |
The local tetrahedral order parameter, $q(z)$, is given by |
345 |
gezelter |
3897 |
\begin{equation} |
346 |
gezelter |
3946 |
q(z) = \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3} |
347 |
|
|
\sum_{j=i+1}^{4} \bigg(\cos\psi_{ikj}+\frac{1}{3}\bigg)^2\Bigg) |
348 |
|
|
\delta(z_{k}-z)\mathrm{d}z \Bigg/ N_z |
349 |
gezelter |
3945 |
\label{eq:qz} |
350 |
gezelter |
3897 |
\end{equation} |
351 |
gezelter |
3945 |
where $\psi_{ikj}$ is the angle formed between the oxygen site on |
352 |
gezelter |
3946 |
central molecule $k$, and the oxygen sites on two of the four closest |
353 |
|
|
molecules, $i$ and $j$. Molecules $i$ and $j$ are further restricted |
354 |
|
|
to lie within the first solvation shell of molecule $k$. $N_z = \int |
355 |
|
|
\delta(z_k - z) \mathrm{d}z$ is a normalization factor to account for |
356 |
|
|
the varying population of molecules within each finite-width bin. The |
357 |
|
|
local tetrahedral order parameter has a range of $(0,1)$, where the |
358 |
|
|
larger values of $q$ indicate a larger degree of tetrahedral ordering |
359 |
|
|
of the local environment. In perfect ice I$_\mathrm{h}$ structures, |
360 |
|
|
the parameter can approach 1 at low temperatures, while in liquid |
361 |
|
|
water, the ordering is significantly less tetrahedral, and values of |
362 |
|
|
$q(z) \approx 0.75$ are more common. |
363 |
plouden |
3909 |
|
364 |
gezelter |
3946 |
To estimate the interfacial width, the system was divided into 100 |
365 |
|
|
bins along the $z$-dimension, and a cutoff radius for the first |
366 |
|
|
solvation shell was set to 3.41 \AA\ . The $q_{z}$ function was |
367 |
|
|
time-averaged to give yield a tetrahedrality profile of the |
368 |
|
|
system. The profile was then fit to a hyperbolic tangent that smoothly |
369 |
|
|
links the liquid and solid states, |
370 |
plouden |
3941 |
\begin{equation}\label{tet_fit} |
371 |
gezelter |
3946 |
q(z) \approx |
372 |
|
|
q_{liq}+\frac{q_{ice}-q_{liq}}{2}\left[\tanh\left(\frac{z-l}{w}\right)-\tanh\left(\frac{z-r}{w}\right)\right]+\beta\left|z- |
373 |
|
|
\frac{r+l}{2}\right|. |
374 |
gezelter |
3897 |
\end{equation} |
375 |
gezelter |
3946 |
Here $q_{liq}$ and $q_{ice}$ are the local tetrahedral order parameter |
376 |
|
|
for the bulk liquid and ice domains, respectively, $w$ is the width of |
377 |
|
|
the interface. $l$ and $r$ are the midpoints of the left and right |
378 |
|
|
interfaces, respectively. The last term in equation \eqref{tet_fit} |
379 |
|
|
accounts for the influence that the weak thermal gradient has on the |
380 |
|
|
tetrahedrality profile in the liquid region. To estimate the |
381 |
gezelter |
3954 |
10\%-90\% widths commonly used in previous studies,\cite{Bryk02} it is |
382 |
|
|
a simple matter to scale the widths obtained from the hyperbolic |
383 |
|
|
tangent fits to obtain $w_{10-90} = 2.1971 \times w$.\cite{Bryk02} |
384 |
gezelter |
3897 |
|
385 |
gezelter |
3946 |
In Figures \ref{fig:bComic} and \ref{fig:pComic} we see the |
386 |
|
|
$z$-coordinate profiles for tetrahedrality, temperature, and the |
387 |
|
|
$x$-component of the velocity for the basal and prismatic interfaces. |
388 |
|
|
The lower panels show the $q(z)$ (black circles) along with the |
389 |
|
|
hyperbolic tangent fits (red lines). In the liquid region, the local |
390 |
|
|
tetrahedral order parameter, $q(z) \approx 0.75$ while in the |
391 |
|
|
crystalline region, $q(z) \approx 0.94$, indicating a more tetrahedral |
392 |
|
|
environment. The vertical dotted lines denote the midpoint of the |
393 |
|
|
interfaces ($r$ and $l$ in equation \eqref{tet_fit}). The weak thermal |
394 |
|
|
gradient applied to the systems in order to keep the interface at |
395 |
gezelter |
3954 |
225$\pm$5K, can be seen in middle panels. The transverse velocity |
396 |
gezelter |
3946 |
profile is shown in the upper panels. It is clear from the upper |
397 |
|
|
panels that water molecules in close proximity to the surface (i.e. |
398 |
|
|
within 10 \AA\ to 15 \AA\ of the interfaces) have transverse |
399 |
|
|
velocities quite close to the velocities within the ice block. There |
400 |
|
|
is no velocity discontinuity at the interface, which indicates that |
401 |
|
|
the shearing of ice/water interfaces occurs in the ``stick'' or |
402 |
|
|
no-slip boundary conditions. |
403 |
gezelter |
3897 |
|
404 |
gezelter |
3914 |
\begin{figure} |
405 |
gezelter |
3946 |
\includegraphics[width=\linewidth]{bComicStrip.pdf} |
406 |
|
|
\caption{\label{fig:bComic} The basal interfaces. Lower panel: the |
407 |
|
|
local tetrahedral order parameter, $q(z)$, (black circles) and the |
408 |
|
|
hyperbolic tangent fit (red line). Middle panel: the imposed |
409 |
|
|
thermal gradient required to maintain a fixed interfacial |
410 |
|
|
temperature. Upper panel: the transverse velocity gradient that |
411 |
|
|
develops in response to an imposed momentum flux. The vertical |
412 |
|
|
dotted lines indicate the locations of the midpoints of the two |
413 |
|
|
interfaces.} |
414 |
gezelter |
3914 |
\end{figure} |
415 |
plouden |
3904 |
|
416 |
gezelter |
3914 |
\begin{figure} |
417 |
gezelter |
3946 |
\includegraphics[width=\linewidth]{pComicStrip.pdf} |
418 |
|
|
\caption{\label{fig:pComic} The prismatic interfaces. Panel |
419 |
|
|
descriptions match those in figure \ref{fig:bComic}} |
420 |
gezelter |
3914 |
\end{figure} |
421 |
|
|
|
422 |
gezelter |
3946 |
From the fits using equation \eqref{tet_fit}, we find the interfacial |
423 |
|
|
width for the basal and prismatic systems to be 3.2$\pm$0.4 \AA\ and |
424 |
|
|
3.6$\pm$0.2 \AA\ , respectively, with no applied momentum flux. Over |
425 |
|
|
the range of shear rates investigated, $0.6 \pm 0.3 \mathrm{ms}^{-1} |
426 |
|
|
\rightarrow 5.3 \pm 0.5 \mathrm{ms}^{-1}$ for the basal system and |
427 |
|
|
$0.9 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 4.5 \pm 0.1 |
428 |
|
|
\mathrm{ms}^{-1}$ for the prismatic, we found no appreciable change in |
429 |
|
|
the interface width. The fit values for the interfacial width ($w$) |
430 |
|
|
over all shear rates contained the values reported above within their |
431 |
|
|
error bars. |
432 |
gezelter |
3914 |
|
433 |
gezelter |
3954 |
\subsubsection{Orientational Dynamics} |
434 |
gezelter |
3946 |
The orientational time correlation function, |
435 |
plouden |
3941 |
\begin{equation}\label{C(t)1} |
436 |
gezelter |
3946 |
C_{2}(t)=\langle P_{2}(\mathbf{u}(0) \cdot \mathbf{u}(t)) \rangle, |
437 |
plouden |
3941 |
\end{equation} |
438 |
gezelter |
3946 |
gives insight into the local dynamic environment around the water |
439 |
|
|
molecules. The rate at which the function decays provides information |
440 |
|
|
about hindered motions and the timescales for relaxation. In |
441 |
|
|
eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial, |
442 |
gezelter |
3948 |
the vector $\mathbf{u}$ is often taken as HOH bisector, although |
443 |
|
|
slightly different behavior can be observed when $\mathbf{u}$ is the |
444 |
gezelter |
3946 |
vector along one of the OH bonds. The angle brackets denote an |
445 |
|
|
ensemble average over all water molecules in a given spatial region. |
446 |
gezelter |
3897 |
|
447 |
gezelter |
3946 |
To investigate the dynamic behavior of water at the ice interfaces, we |
448 |
|
|
have computed $C_{2}(z,t)$ for molecules that are present within a |
449 |
|
|
particular slab along the $z$- axis at the initial time. The change |
450 |
|
|
in the decay behavior as a function of the $z$ coordinate is another |
451 |
gezelter |
3948 |
measure of the change of how the local environment changes across the |
452 |
|
|
ice/water interface. To compute these correlation functions, each of |
453 |
|
|
the 0.5 ns simulations was followed by a shorter 200 ps microcanonical |
454 |
|
|
(NVE) simulation in which the positions and orientations of every |
455 |
|
|
molecule in the system were recorded every 0.1 ps. The systems were |
456 |
gezelter |
3954 |
then divided into 30 bins along the $z$-axis and $C_2(t)$ was |
457 |
|
|
evaluated for each bin. |
458 |
plouden |
3919 |
|
459 |
gezelter |
3948 |
In simulations of water at biological interfaces, Furse {\em et al.} |
460 |
|
|
fit $C_2(t)$ functions for water with triexponential |
461 |
|
|
functions,\cite{Furse08} where the three components of the decay |
462 |
|
|
correspond to a fast (<200 fs) reorientational piece driven by the |
463 |
|
|
restoring forces of existing hydrogen bonds, a middle (on the order of |
464 |
|
|
several ps) piece describing the large angle jumps that occur during |
465 |
|
|
the breaking and formation of new hydrogen bonds,and a slow (on the |
466 |
|
|
order of tens of ps) contribution describing the translational motion |
467 |
|
|
of the molecules. The model for orientational decay presented |
468 |
|
|
recently by Laage and Hynes {\em et al.}\cite{Laage08,Laage11} also |
469 |
|
|
includes three similar decay constants, although two of the time |
470 |
|
|
constants are linked, and the resulting decay curve has two parameters |
471 |
|
|
governing the dynamics of decay. |
472 |
|
|
|
473 |
|
|
In our ice/water interfaces, we are at substantially lower |
474 |
|
|
temperatures, and the water molecules are further perturbed by the |
475 |
|
|
presence of the ice phase nearby. We have obtained the most |
476 |
|
|
reasonable fits using triexponential functions with three distinct |
477 |
gezelter |
3954 |
time domains, as well as a constant piece to account for the water |
478 |
gezelter |
3948 |
stuck in the ice phase that does not experience any long-time |
479 |
|
|
orientational decay, |
480 |
gezelter |
3897 |
\begin{equation} |
481 |
gezelter |
3946 |
C_{2}(t) \approx a e^{-t/\tau_\mathrm{short}} + b e^{-t/\tau_\mathrm{middle}} + c |
482 |
|
|
e^{-t/\tau_\mathrm{long}} + (1-a-b-c) |
483 |
gezelter |
3897 |
\end{equation} |
484 |
gezelter |
3948 |
Average values for the three decay constants (and error estimates) |
485 |
|
|
were obtained for each bin. In figures \ref{fig:basal_Tau_comic_strip} |
486 |
|
|
and \ref{fig:prismatic_Tau_comic_strip}, the three orientational decay |
487 |
|
|
times are shown as a function of distance from the center of the ice |
488 |
|
|
slab. |
489 |
gezelter |
3897 |
|
490 |
plouden |
3941 |
\begin{figure} |
491 |
gezelter |
3946 |
\includegraphics[width=\linewidth]{basal_Tau_comic_strip.pdf} |
492 |
gezelter |
3949 |
\caption{\label{fig:basal_Tau_comic_strip} The three decay constants |
493 |
|
|
of the orientational time correlation function, $C_2(t)$, for water |
494 |
|
|
as a function of distance from the center of the ice slab. The |
495 |
|
|
dashed line indicates the location of the basal face (as determined |
496 |
|
|
from the tetrahedrality order parameter). The moderate and long |
497 |
|
|
time contributions slow down close to the interface which would be |
498 |
gezelter |
3954 |
expected under reorganizations that involve large motions of the |
499 |
gezelter |
3949 |
molecules (e.g. frame-reorientations and jumps). The observed |
500 |
|
|
speed-up in the short time contribution is surprising, but appears |
501 |
|
|
to reflect the restricted motion of librations closer to the |
502 |
|
|
interface.} |
503 |
plouden |
3941 |
\end{figure} |
504 |
gezelter |
3897 |
|
505 |
gezelter |
3914 |
\begin{figure} |
506 |
gezelter |
3946 |
\includegraphics[width=\linewidth]{prismatic_Tau_comic_strip.pdf} |
507 |
gezelter |
3949 |
\caption{\label{fig:prismatic_Tau_comic_strip} |
508 |
|
|
Decay constants for $C_2(t)$ at the prismatic interface. Panel |
509 |
|
|
descriptions match those in figure \ref{fig:basal_Tau_comic_strip}.} |
510 |
gezelter |
3914 |
\end{figure} |
511 |
plouden |
3904 |
|
512 |
gezelter |
3949 |
Figures \ref{fig:basal_Tau_comic_strip} and |
513 |
|
|
\ref{fig:prismatic_Tau_comic_strip} show the three decay constants for |
514 |
|
|
the orientational time correlation function for water at varying |
515 |
|
|
displacements from the center of the ice slab for both the basal and |
516 |
|
|
prismatic interfaces. The vertical dotted lines indicate the |
517 |
|
|
locations of the midpoints of the interfaces as determined by the |
518 |
|
|
tetrahedrality fits. In the liquid regions, $\tau_{middle}$ and |
519 |
|
|
$\tau_{long}$ have consistent values around 3-4 ps and 20-40 ps, |
520 |
|
|
respectively, and increase in value approaching the interface. |
521 |
|
|
According to the jump model of Laage and Hynes {\em et |
522 |
|
|
al.},\cite{Laage08,Laage11} $\tau_{middle}$ corresponds to the |
523 |
|
|
breaking and making of hydrogen bonds and $\tau_{long}$ is explained |
524 |
|
|
with translational motion of the molecules (i.e. frame reorientation). |
525 |
|
|
The shortest of the three decay constants, the librational time |
526 |
|
|
$\tau_\mathrm{short}$ has a value of about 70 fs in the liquid region, |
527 |
|
|
and decreases in value approaching the interface. The observed |
528 |
|
|
speed-up in the short time contribution is surprising, but appears to |
529 |
|
|
reflect the restricted motion of librations closer to the interface. |
530 |
plouden |
3937 |
|
531 |
gezelter |
3949 |
The control systems (with no applied momentum flux) are shown with |
532 |
|
|
black symbols in figs. \ref{fig:basal_Tau_comic_strip} and |
533 |
|
|
\ref{fig:prismatic_Tau_comic_strip}, while those obtained while a |
534 |
|
|
shear was active are shown in red. |
535 |
plouden |
3937 |
|
536 |
gezelter |
3952 |
Two notable features deserve clarification. First, there are |
537 |
|
|
nearly-constant liquid-state values for $\tau_{short}$, |
538 |
gezelter |
3949 |
$\tau_{middle}$, and $\tau_{long}$ at large displacements from the |
539 |
gezelter |
3952 |
interface. Second, there appears to be a single distance, $d_{basal}$ |
540 |
|
|
or $d_{prismatic}$, from the interface at which all three decay times |
541 |
|
|
begin to deviate from their bulk liquid values. We find these |
542 |
|
|
distances to be approximately 15 \AA\ and 8 \AA\ , respectively, |
543 |
|
|
although significantly finer binning of the $C_2(t)$ data would be |
544 |
|
|
necessary to provide better estimates of a ``dynamic'' interfacial |
545 |
|
|
thickness. |
546 |
plouden |
3937 |
|
547 |
gezelter |
3952 |
Beaglehole and Wilson have measured the ice/water interface using |
548 |
|
|
ellipsometry and find a thickness of approximately 10 \AA\ for both |
549 |
|
|
the basal and prismatic faces.\cite{Beaglehole93} Structurally, we |
550 |
|
|
have found the basal and prismatic interfacial width to be 3.2$\pm$0.4 |
551 |
|
|
\AA\ and 3.6$\pm$0.2 \AA\ . However, decomposition of the spatial |
552 |
|
|
dependence of the decay times of $C_2(t)$ indicates that a somewhat |
553 |
|
|
thicker interfacial region exists in which the orientational dynamics |
554 |
|
|
of the water molecules begin to resemble the trapped interfacial water |
555 |
|
|
more than the surrounding liquid. |
556 |
plouden |
3937 |
|
557 |
gezelter |
3952 |
Our results indicate that the dynamics of the water molecules within |
558 |
|
|
$d_{basal}$ and $d_{prismatic}$ are being significantly perturbed by |
559 |
|
|
the interface, even though the structural width of the interface via |
560 |
|
|
analysis of the tetrahedrality profile indicates that bulk liquid |
561 |
|
|
structure of water is recovered after about 4 \AA\ from the edge of |
562 |
|
|
the ice. |
563 |
|
|
|
564 |
plouden |
3941 |
\subsection{Coefficient of Friction of the Interface} |
565 |
gezelter |
3954 |
As liquid water flows over an ice interface, there is a distance from |
566 |
|
|
the structural interface where bulk-like hydrodynamics are recovered. |
567 |
|
|
Bocquet and Barrat constructed a theory for the hydrodynamic boundary |
568 |
|
|
parameters, which include the slipping length |
569 |
gezelter |
3952 |
$\left(\delta_\mathrm{wall}\right)$ of this boundary layer and the |
570 |
|
|
``hydrodynamic position'' of the boundary |
571 |
gezelter |
3954 |
$\left(z_\mathrm{wall}\right)$.\cite{PhysRevLett.70.2726,PhysRevE.49.3079} |
572 |
|
|
This last parameter is the location (relative to a solid surface) |
573 |
gezelter |
3952 |
where the bulk-like behavior is recovered. Work by Mundy {\it et al.} |
574 |
gezelter |
3954 |
has helped to combine these parameters into a liquid-solid friction |
575 |
gezelter |
3952 |
coefficient, which quantifies the resistance to pulling the solid |
576 |
gezelter |
3954 |
interface through a liquid,\cite{Mundy1997305} |
577 |
gezelter |
3952 |
\begin{equation} |
578 |
|
|
\lambda_\mathrm{wall} = \frac{\eta}{\delta_\mathrm{wall}}. |
579 |
|
|
\end{equation} |
580 |
gezelter |
3954 |
This expression is nearly identical to one provided by Pit {\it et |
581 |
|
|
al.} for the solid-liquid friction of an interface,\cite{Pit99} |
582 |
plouden |
3942 |
\begin{equation}\label{Pit} |
583 |
gezelter |
3952 |
\lambda=\frac{\eta}{\delta} |
584 |
plouden |
3942 |
\end{equation} |
585 |
gezelter |
3952 |
where $\delta$ is the slip length for the liquid measured at the |
586 |
|
|
location of the interface itself. |
587 |
|
|
|
588 |
|
|
In both of these expressions, $\eta$ is the shear viscosity of the |
589 |
|
|
bulk-like region of the liquid, a quantity which is easily obtained in |
590 |
|
|
VSS-RNEMD simulations by fitting the velocity profile in the region |
591 |
gezelter |
3954 |
far from the surface.\cite{Kuang12} Assuming linear response in the |
592 |
gezelter |
3952 |
bulk-like region, |
593 |
plouden |
3942 |
\begin{equation}\label{Kuang} |
594 |
gezelter |
3952 |
j_{z}(p_{x})=-\eta \left(\frac{\partial v_{x}}{\partial z}\right) |
595 |
plouden |
3942 |
\end{equation} |
596 |
gezelter |
3952 |
Substituting this result into eq. \eqref{Pit}, we can estimate the |
597 |
|
|
solid-liquid coefficient using the slip length, |
598 |
plouden |
3941 |
\begin{equation} |
599 |
gezelter |
3954 |
\lambda=-\frac{j_{z}(p_{x})} {\left(\frac{\partial v_{x}}{\partial |
600 |
|
|
z}\right) \delta} |
601 |
plouden |
3941 |
\end{equation} |
602 |
plouden |
3937 |
|
603 |
gezelter |
3954 |
For ice / water interfaces, the boundary conditions are markedly |
604 |
|
|
no-slip, so projecting the bulk liquid state velocity profile yields a |
605 |
|
|
negative slip length. This length is the difference between the |
606 |
|
|
structural edge of the ice (determined by the tetrahedrality profile) |
607 |
|
|
and the location where the projected velocity of the bulk liquid |
608 |
|
|
intersects the solid phase velocity (see Figure |
609 |
|
|
\ref{fig:delta_example}). The coefficients of friction for the basal |
610 |
|
|
and the prismatic facets are found to be (0.07$\pm$0.01) amu |
611 |
gezelter |
3952 |
fs\textsuperscript{-1} and (0.032$\pm$0.007) amu |
612 |
gezelter |
3954 |
fs\textsuperscript{-1}, respectively. These results may seem |
613 |
|
|
surprising as the basal face is smoother than the prismatic with only |
614 |
|
|
small undulations of the oxygen positions, while the prismatic surface |
615 |
|
|
has deep corrugated channels. The applied momentum flux used in our |
616 |
|
|
simulations is parallel to these channels, however, and this results |
617 |
|
|
in a flow of water in the same direction as the corrugations, allowing |
618 |
|
|
water molecules to pass through the channels during the shear. |
619 |
plouden |
3939 |
|
620 |
plouden |
3942 |
\begin{figure} |
621 |
gezelter |
3946 |
\includegraphics[width=\linewidth]{delta_example.pdf} |
622 |
gezelter |
3952 |
\caption{\label{fig:delta_example} Determining the (negative) slip |
623 |
|
|
length ($\delta$) for the ice-water interfaces (which have decidedly |
624 |
|
|
non-slip behavior). This length is the difference between the |
625 |
|
|
structural edge of the ice (determined by the tetrahedrality |
626 |
|
|
profile) and the location where the projected velocity of the bulk |
627 |
|
|
liquid (dashed red line) intersects the solid phase velocity (solid |
628 |
|
|
black line). The dotted line indicates the location of the ice as |
629 |
|
|
determined by the tetrahedrality profile.} |
630 |
plouden |
3942 |
\end{figure} |
631 |
|
|
|
632 |
|
|
|
633 |
gezelter |
3897 |
\section{Conclusion} |
634 |
gezelter |
3952 |
We have simulated the basal and prismatic facets of an SPC/E model of |
635 |
|
|
the ice I$_\mathrm{h}$ / water interface. Using VSS-RNEMD, the ice |
636 |
|
|
was sheared relative to the liquid while simultaneously being exposed |
637 |
|
|
to a weak thermal gradient which kept the interface at a stable |
638 |
|
|
temperature. Calculation of the local tetrahedrality order parameter |
639 |
|
|
has shown an apparent independence of the interfacial width on the |
640 |
|
|
shear rate. This width was found to be 3.2$\pm$0.4 \AA\ and |
641 |
|
|
3.6$\pm$0.2 \AA\ for the basal and prismatic systems, respectively. |
642 |
gezelter |
3897 |
|
643 |
gezelter |
3952 |
Orientational time correlation functions were calculated at varying |
644 |
|
|
displacements from the interface, and were found to be similarly |
645 |
|
|
independent of the applied momentum flux. The short decay due to the |
646 |
|
|
restoring forces of existing hydrogen bonds decreased close to the |
647 |
|
|
interface, while the longer-time decay constants increased in close |
648 |
|
|
proximity to the interface. There is also an apparent dynamic |
649 |
|
|
interface width, $d_{basal}$ and $d_{prismatic}$, at which these |
650 |
gezelter |
3954 |
deviations from bulk liquid values begin. We found $d_{basal}$ and |
651 |
|
|
$d_{prismatic}$ to be approximately 15 \AA\ and 8 \AA\ . This implies |
652 |
|
|
that the dynamics of water molecules which have similar structural |
653 |
|
|
environments to liquid phase molecules are dynamically perturbed by |
654 |
|
|
the presence of the ice interface. |
655 |
gezelter |
3952 |
|
656 |
gezelter |
3954 |
The coefficient of liquid-solid friction for each of the facets was |
657 |
|
|
also determined. They were found to be (0.07$\pm$0.01) amu |
658 |
gezelter |
3952 |
fs\textsuperscript{-1} and (0.032$\pm$0.007) amu |
659 |
|
|
fs\textsuperscript{-1} for the basal and prismatic facets |
660 |
|
|
respectively. We attribute the large difference between the two |
661 |
|
|
friction coefficients to the direction of the shear and to the surface |
662 |
|
|
structure of the crystal facets. |
663 |
|
|
|
664 |
gezelter |
3914 |
\begin{acknowledgement} |
665 |
|
|
Support for this project was provided by the National Science |
666 |
|
|
Foundation under grant CHE-0848243. Computational time was provided |
667 |
|
|
by the Center for Research Computing (CRC) at the University of |
668 |
|
|
Notre Dame. |
669 |
|
|
\end{acknowledgement} |
670 |
gezelter |
3897 |
|
671 |
gezelter |
3914 |
\newpage |
672 |
|
|
\bibstyle{achemso} |
673 |
plouden |
3909 |
\bibliography{iceWater} |
674 |
gezelter |
3897 |
|
675 |
gezelter |
3914 |
\begin{tocentry} |
676 |
|
|
\begin{wrapfigure}{l}{0.5\textwidth} |
677 |
|
|
\begin{center} |
678 |
|
|
\includegraphics[width=\linewidth]{SystemImage.png} |
679 |
|
|
\end{center} |
680 |
|
|
\end{wrapfigure} |
681 |
gezelter |
3954 |
The cell used to simulate liquid-solid shear in ice I$_\mathrm{h}$ / |
682 |
|
|
water interfaces. Velocity gradients were applied using the velocity |
683 |
|
|
shearing and scaling variant of reverse non-equilibrium molecular |
684 |
|
|
dynamics (VSS-RNEMD) with a weak thermal gradient to prevent melting. |
685 |
|
|
The interface width is relatively robust in both structual and dynamic |
686 |
|
|
measures as a function of the applied shear. |
687 |
gezelter |
3914 |
\end{tocentry} |
688 |
gezelter |
3897 |
|
689 |
gezelter |
3914 |
\end{document} |
690 |
gezelter |
3897 |
|
691 |
plouden |
3924 |
%************************************************************** |
692 |
|
|
%Non-mass weighted slopes (amu*Angstroms^-2 * fs^-1) |
693 |
|
|
% basal: slope=0.090677616, error in slope = 0.003691743 |
694 |
|
|
%prismatic: slope = 0.050101506, error in slope = 0.001348181 |
695 |
|
|
%Mass weighted slopes (Angstroms^-2 * fs^-1) |
696 |
|
|
%basal slope = 4.76598E-06, error in slope = 1.94037E-07 |
697 |
|
|
%prismatic slope = 3.23131E-06, error in slope = 8.69514E-08 |
698 |
|
|
%************************************************************** |