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 |
|
|
|
14 |
gezelter |
3945 |
\title{Solid-liquid friction at ice-I$_\mathrm{h}$ / 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 |
|
|
basal and prismatic facets of an ice I$_\mathrm{h}$ / water |
30 |
|
|
interface when the solid phase is being drawn through the liquid |
31 |
|
|
(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 |
|
|
width was found to be independent of relative velocity of the ice |
37 |
|
|
and liquid layers over a wide range of shear rates. Decays of |
38 |
|
|
molecular orientational time correlation functions for gave very |
39 |
|
|
similar estimates for the width of the interfaces, although the |
40 |
|
|
short- and longer-time decay components of the orientational |
41 |
|
|
correlation functions behave differently closer to the interface. |
42 |
|
|
Although both facets of ice are in ``stick'' boundary conditions in |
43 |
|
|
liquid water, the solid-liquid friction coefficient was found to be |
44 |
|
|
different for the basal and prismatic 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 |
|
|
interface for the SPC, SPC/E, CF1, and TIP4P models for |
71 |
|
|
water.\cite{Karim87,Karim88,Karim90,Hayward01,Bryk02,Hayward02,Gay02} |
72 |
|
|
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 |
3945 |
The structure of ice I$_\mathrm{h}$ is very 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 |
|
|
the prismatic, $\{2~0~\bar{2}~1\}$, faces. Ice I$_\mathrm{h}$ is |
120 |
|
|
normally proton disordered in bulk crystals, although the surfaces |
121 |
|
|
probably have a preference for proton ordering along strips of |
122 |
|
|
dangling H-atoms and Oxygen lone pairs.\cite{Buch:2008fk} |
123 |
gezelter |
3914 |
|
124 |
|
|
For small simulated ice interfaces, it is useful to have a |
125 |
|
|
proton-ordered, but zero-dipole crystal that exposes these strips of |
126 |
gezelter |
3945 |
dangling H-atoms and lone pairs. When placing another material in |
127 |
|
|
contact with one of the ice crystalline planes, it is useful to have |
128 |
|
|
an orthorhombic (rectangular) box. A recent paper by Hirsch and |
129 |
|
|
Ojam\"{a}e describes how to create proton-ordered bulk ice |
130 |
|
|
I$_\mathrm{h}$ in alternative orthorhombic cells.\cite{Hirsch04} |
131 |
gezelter |
3914 |
|
132 |
plouden |
3941 |
We are using Hirsch and Ojam\"{a}e's structure 6 which is an |
133 |
gezelter |
3914 |
orthorhombic cell ($P2_12_12_1$) that produces a proton-ordered |
134 |
gezelter |
3945 |
version of ice I$_\mathrm{h}$. Table \ref{tab:equiv} contains a mapping between |
135 |
gezelter |
3914 |
the Miller indices in the P$6_3/mmc$ crystal system and those in the |
136 |
|
|
Hirsch and Ojam\"{a}e $P2_12_12_1$ system. |
137 |
|
|
|
138 |
|
|
\begin{wraptable}{r}{3.5in} |
139 |
gezelter |
3945 |
\caption{Mapping between the Miller indices of four facets of ice in |
140 |
|
|
the $P6_3/mmc$ crystal system to the orthorhombic $P2_12_12_1$ |
141 |
|
|
system in reference \protect\cite{Hirsch04}} |
142 |
|
|
\label{tab:equiv} |
143 |
gezelter |
3914 |
\begin{tabular}{|ccc|} \hline |
144 |
|
|
& hexagonal & orthorhombic \\ |
145 |
|
|
& ($P6_3/mmc$) & ($P2_12_12_1$) \\ |
146 |
|
|
crystal face & Miller indices & equivalent \\ \hline |
147 |
|
|
basal & $\{0~0~0~1\}$ & $\{0~0~1\}$ \\ |
148 |
|
|
prism & $\{1~0~\bar{1}~0\}$ & $\{1~0~0\}$ \\ |
149 |
|
|
secondary prism & $\{1~1~\bar{2}~0\}$ & $\{1~3~0\}$ \\ |
150 |
|
|
pryamidal & $\{2~0~\bar{2}~1\}$ & $\{2~0~1\}$ \\ \hline |
151 |
|
|
\end{tabular} |
152 |
|
|
\end{wraptable} |
153 |
|
|
|
154 |
|
|
Structure 6 from the Hirsch and Ojam\"{a}e paper has lattice |
155 |
|
|
parameters $a = 4.49225$, $b = 7.78080$, $c = 7.33581$ and two water |
156 |
gezelter |
3945 |
molecules whose atoms reside at fractional coordinates given in table |
157 |
|
|
\ref{tab:p212121}. To construct the basal and prismatic interfaces, |
158 |
|
|
these crystallographic coordinates were used to construct an |
159 |
|
|
orthorhombic unit cell which was then replicated in all three |
160 |
|
|
dimensions yielding a proton-ordered block of ice I$_{h}$. To expose |
161 |
|
|
the desired face, the orthorhombic representation was then cut along |
162 |
|
|
the ($001$) or ($100$) planes for the basal and prismatic faces |
163 |
|
|
respectively. The resulting block was rotated so that the exposed |
164 |
|
|
faces were aligned with the $z$-dimension normal to the exposed face. |
165 |
|
|
The block was then cut along two perpendicular directions in a way |
166 |
|
|
that allowed for perfect periodic replication in the $x$ and $y$ axes, |
167 |
|
|
creating a slab with either the basal or prismatic faces exposed along |
168 |
|
|
the $z$ axis. The slab was then replicated in the $x$ and $y$ |
169 |
|
|
dimensions until a desired sample size was obtained. |
170 |
gezelter |
3914 |
|
171 |
|
|
\begin{wraptable}{r}{3.25in} |
172 |
gezelter |
3945 |
\caption{Fractional coordinates for water in the orthorhombic |
173 |
|
|
$P2_12_12_1$ system for ice I$_\mathrm{h}$ in reference |
174 |
|
|
\protect\cite{Hirsch04}} |
175 |
|
|
\label{tab:p212121} |
176 |
|
|
\begin{tabular}{|cccc|} \hline |
177 |
|
|
atom type & x & y & z \\ \hline |
178 |
|
|
O & 0.75 & 0.1667 & 0.4375 \\ |
179 |
|
|
H & 0.5735 & 0.2202 & 0.4836 \\ |
180 |
|
|
H & 0.7420 & 0.0517 & 0.4836 \\ |
181 |
|
|
O & 0.25 & 0.6667 & 0.4375 \\ |
182 |
|
|
H & 0.2580 & 0.6693 & 0.3071 \\ |
183 |
|
|
H & 0.4265 & 0.7255 & 0.4756 \\ \hline |
184 |
gezelter |
3914 |
\end{tabular} |
185 |
|
|
\end{wraptable} |
186 |
|
|
|
187 |
|
|
Although experimental solid/liquid coexistant temperature under normal |
188 |
|
|
pressure are close to 273K, Haymet \emph{et al.} have done extensive |
189 |
|
|
work on characterizing the ice/water |
190 |
|
|
interface.\cite{Karim88,Karim90,Hayward01,Bryk02,Hayward02} They have |
191 |
gezelter |
3945 |
found that for the SPC/E water model,\cite{Berendsen87} which is also |
192 |
|
|
used in this study, the ice/water interface is most stable at |
193 |
|
|
225$\pm$5K.\cite{Bryk02} Therefore, we created our ice / water |
194 |
|
|
interfaces, utilizing a box of liquid water that had the same |
195 |
|
|
dimensions in $x$ and $y$ was equilibrated at 225 K and 1 atm of |
196 |
|
|
pressure in the NPAT ensemble (with the $z$ axis allowed to fluctuate |
197 |
|
|
to equilibrate to the correct pressure). The liquid and solid systems |
198 |
|
|
were combined by carving out any water molecule from the liquid |
199 |
|
|
simulation cell that was within 3 \AA\ of any atom in the ice slab. |
200 |
gezelter |
3914 |
|
201 |
|
|
Molecular translation and orientational restraints were applied in the |
202 |
|
|
early stages of equilibration to prevent melting of the ice slab. |
203 |
|
|
These restraints were removed during NVT equilibration, well before |
204 |
gezelter |
3945 |
data collection was carried out. |
205 |
gezelter |
3914 |
|
206 |
gezelter |
3918 |
\subsection{Shearing ice / water interfaces without bulk melting} |
207 |
|
|
|
208 |
gezelter |
3945 |
As a solid is dragged through a liquid, there is frictional heating |
209 |
|
|
that will act to melt the interface. To study the behavior of the |
210 |
|
|
interface under a shear stress without causing the interface to melt, |
211 |
|
|
it is necessary to apply a weak thermal gradient along with the |
212 |
|
|
momentum gradient. This can be accomplished using he velocity |
213 |
|
|
shearing and scaling (VSS) variant of reverse non-equilibrium |
214 |
|
|
molecular dynamics (RNEMD), which utilizes a series of simultaneous |
215 |
|
|
velocity exchanges between two regions within the simulation |
216 |
|
|
cell.\cite{Kuang12} One of these regions is centered within the ice |
217 |
|
|
slab, while the other is centrally located in the liquid phase |
218 |
gezelter |
3918 |
region. VSS-RNEMD provides a set of conservation constraints for |
219 |
|
|
simultaneously creating either a momentum flux or a thermal flux (or |
220 |
|
|
both) between the two slabs. Satisfying the constraint equations |
221 |
|
|
ensures that the new configurations are sampled from the same NVE |
222 |
gezelter |
3945 |
ensemble as before the VSS move. |
223 |
gezelter |
3918 |
|
224 |
|
|
The VSS moves are applied periodically to scale and shift the particle |
225 |
|
|
velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and |
226 |
|
|
$C$) which are separated by half of the simulation box, |
227 |
|
|
\begin{displaymath} |
228 |
|
|
\begin{array}{rclcl} |
229 |
|
|
|
230 |
|
|
& \underline{\mathrm{shearing}} & & |
231 |
|
|
\underline{~~~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~~~} \\ |
232 |
|
|
\mathbf{v}_i \leftarrow & |
233 |
|
|
\mathbf{a}_c\ & + & c\cdot\left(\mathbf{v}_i - \langle\mathbf{v}_c |
234 |
|
|
\rangle\right) + \langle\mathbf{v}_c\rangle \\ |
235 |
|
|
\mathbf{v}_j \leftarrow & |
236 |
|
|
\mathbf{a}_h & + & h\cdot\left(\mathbf{v}_j - \langle\mathbf{v}_h |
237 |
|
|
\rangle\right) + \langle\mathbf{v}_h\rangle . |
238 |
|
|
|
239 |
|
|
\end{array} |
240 |
|
|
\end{displaymath} |
241 |
|
|
Here $\langle\mathbf{v}_c\rangle$ and $\langle\mathbf{v}_h\rangle$ are |
242 |
|
|
the center of mass velocities in the $C$ and $H$ slabs, respectively. |
243 |
|
|
Within the two slabs, particles receive incremental changes or a |
244 |
|
|
``shear'' to their velocities. The amount of shear is governed by the |
245 |
|
|
imposed momentum flux, $\mathbf{j}_z(\mathbf{p})$ |
246 |
|
|
\begin{eqnarray} |
247 |
|
|
\mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\ |
248 |
|
|
\mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2} |
249 |
|
|
\end{eqnarray} |
250 |
|
|
where $M_{\{c,h\}}$ is the total mass of particles within each of the |
251 |
|
|
slabs and $\Delta t$ is the interval between two separate operations. |
252 |
|
|
|
253 |
|
|
To simultaneously impose a thermal flux ($J_z$) between the slabs we |
254 |
|
|
use energy conservation constraints, |
255 |
|
|
\begin{eqnarray} |
256 |
|
|
K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\mathbf{v}_c |
257 |
|
|
\rangle^2) + \frac{1}{2}M_c (\langle \mathbf{v}_c \rangle + \mathbf{a}_c)^2 \label{vss3}\\ |
258 |
|
|
K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\mathbf{v}_h |
259 |
|
|
\rangle^2) + \frac{1}{2}M_h (\langle \mathbf{v}_h \rangle + |
260 |
|
|
\mathbf{a}_h)^2 \label{vss4}. |
261 |
|
|
\label{constraint} |
262 |
|
|
\end{eqnarray} |
263 |
|
|
Simultaneous solution of these quadratic formulae for the scaling |
264 |
|
|
coefficients, $c$ and $h$, will ensure that the simulation samples from |
265 |
|
|
the original microcanonical (NVE) ensemble. Here $K_{\{c,h\}}$ is the |
266 |
|
|
instantaneous translational kinetic energy of each slab. At each time |
267 |
|
|
interval, it is a simple matter to solve for $c$, $h$, $\mathbf{a}_c$, |
268 |
|
|
and $\mathbf{a}_h$, subject to the imposed momentum flux, |
269 |
|
|
$j_z(\mathbf{p})$, and thermal flux, $J_z$, values. Since the VSS |
270 |
|
|
operations do not change the kinetic energy due to orientational |
271 |
|
|
degrees of freedom or the potential energy of a system, configurations |
272 |
|
|
after the VSS move have exactly the same energy (and linear |
273 |
|
|
momentum) as before the move. |
274 |
|
|
|
275 |
|
|
As the simulation progresses, the VSS moves are performed on a regular |
276 |
|
|
basis, and the system develops a thermal and/or velocity gradient in |
277 |
gezelter |
3945 |
response to the applied flux. In a bulk material, it is quite simple |
278 |
gezelter |
3918 |
to use the slope of the temperature or velocity gradients to obtain |
279 |
gezelter |
3945 |
either the thermal conductivity or shear viscosity. |
280 |
gezelter |
3918 |
|
281 |
|
|
The VSS-RNEMD approach is versatile in that it may be used to |
282 |
|
|
implement thermal and shear transport simultaneously. Perturbations |
283 |
|
|
of velocities away from the ideal Maxwell-Boltzmann distributions are |
284 |
|
|
minimal, as is thermal anisotropy. This ability to generate |
285 |
|
|
simultaneous thermal and shear fluxes has been previously utilized to |
286 |
|
|
map out the shear viscosity of SPC/E water over a wide range of |
287 |
|
|
temperatures (90~K) with a single 1 ns simulation.\cite{Kuang12} |
288 |
|
|
|
289 |
gezelter |
3945 |
For this work, we are using the VSS-RNEMD method primarily to generate |
290 |
|
|
a shear between the ice slab and the liquid phase, while using a weak |
291 |
|
|
thermal gradient to maintaining the interface at the 225K target |
292 |
|
|
value. This will insure minimal melting of the bulk ice phase and |
293 |
|
|
allows us to control the exact temperature of the interface. |
294 |
gezelter |
3918 |
|
295 |
gezelter |
3897 |
\subsection{Computational Details} |
296 |
gezelter |
3945 |
All simulations were performed using OpenMD,\cite{OOPSE,openmd} with a |
297 |
|
|
time step of 2 fs and periodic boundary conditions in all three |
298 |
|
|
dimensions. Electrostatics were handled using the damped-shifted |
299 |
|
|
force real-space electrostatic kernel.\cite{Ewald} The systems were |
300 |
|
|
divided into 100 bins along the $z$-axis for the VSS-RNEMD moves, |
301 |
|
|
which were attempted every 50 fs. |
302 |
gezelter |
3897 |
|
303 |
gezelter |
3945 |
The interfaces were equilibrated for a total of 10 ns at equilibrium |
304 |
|
|
conditions before being exposed to either a shear or thermal gradient. |
305 |
|
|
This consisted of 5 ns under a constant temperature (NVT) integrator |
306 |
|
|
set to 225K followed by 5 ns under a microcanonical integrator. Weak |
307 |
|
|
thermal gradients were allowed to develop using the VSS-RNEMD (NVE) |
308 |
|
|
integrator using a a small thermal flux ($-2.0\times 10^{-6}$ |
309 |
|
|
kcal/mol/\AA$^2$/fs) for a duration of 5 ns to allow the gradient to |
310 |
|
|
stabilize. The resulting temperature gradient was less than 5K over |
311 |
|
|
the entire 1 nm box length, which was sufficient to keep the |
312 |
|
|
temperature at the interface within $\pm 1$ K of the 225K target. |
313 |
gezelter |
3897 |
|
314 |
gezelter |
3945 |
Velocity gradients were then imposed using the VSS-RNEMD (NVE) |
315 |
|
|
integrator with a range of momentum fluxes. These gradients were |
316 |
|
|
allowed to stabilize for 1 ns before data collection began. Once |
317 |
|
|
established, four successive 0.5 ns runs were performed for each shear |
318 |
|
|
rate. During these simulations, snapshots of the system were taken |
319 |
|
|
every 1 ps, and statistics on the structure and dynamics in each bin |
320 |
|
|
were accumulated throughout the simulations. |
321 |
|
|
|
322 |
plouden |
3941 |
\section{Results and discussion} |
323 |
|
|
|
324 |
gezelter |
3897 |
\subsection{Measuring the Width of the Interface} |
325 |
gezelter |
3945 |
Any order parameter or correlation function that varies across the |
326 |
|
|
interface from a bulk liquid to a solid can be used as a measure of |
327 |
|
|
the width of the interface. However, because VSS-RNEMD imposes a |
328 |
|
|
lateral flow, parameters that depend on translational motion of the |
329 |
|
|
molecules (e.g. the diffusion constant) may be artifically skewed by |
330 |
|
|
the RNEMD moves. A structural parameter like a radial distribution |
331 |
|
|
function is not influenced by the RNEMD perturbations to the same |
332 |
|
|
degree. Here, we have used the local tetraherdal order parameter as |
333 |
|
|
described by Kumar\cite{Kumar09} and Errington\cite{Errington01} as a |
334 |
|
|
measure of the interfacial width. |
335 |
plouden |
3936 |
|
336 |
gezelter |
3945 |
The local tetrahedral order parameter, $q(z)$, is given by |
337 |
gezelter |
3897 |
\begin{equation} |
338 |
gezelter |
3945 |
q(z) \equiv \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3} |
339 |
|
|
\sum_{j=i+1}^{4} \bigg[\cos\psi_{ikj}+\frac{1}{3}\bigg]^2\Bigg) |
340 |
|
|
\delta(z_{k}-z)\mathrm{d}z |
341 |
|
|
\label{eq:qz} |
342 |
gezelter |
3897 |
\end{equation} |
343 |
gezelter |
3945 |
where $\psi_{ikj}$ is the angle formed between the oxygen site on |
344 |
|
|
molecule $k$, and the oxygen sites on its two closest neighbors, |
345 |
|
|
molecules $i$ and $j$. The local tetrahedral order parameter function |
346 |
|
|
has a range of (0,1), where the larger the value $q$ has the more |
347 |
|
|
tetrahedral the ordering of the local environment is. A $q$ value of |
348 |
|
|
one describes a perfectly tetrahedral environment relative to it and |
349 |
|
|
its four nearest neighbors, and the parameter's value decreases as the |
350 |
|
|
local ordering becomes less tetrahedral. Equation \ref{eq:qz} |
351 |
|
|
describes a $z$-binned tetrahedral order parameter in which the $z$ |
352 |
|
|
coordinate of the central molecule is used to give a spatial |
353 |
|
|
description of the local orientational ordering. |
354 |
plouden |
3909 |
|
355 |
gezelter |
3945 |
The system was divided into 100 bins along the $z$-dimension, and a |
356 |
|
|
cutoff radius for the neighboring molecules was set to 3.41 \AA\ . |
357 |
|
|
The $q_{z}$ values for each snapshot were then averaged to give a |
358 |
|
|
tetrahedrality profile of the system about the $z$-dimension. The |
359 |
|
|
profile was then fit with a hyperbolic tangent function given by |
360 |
plouden |
3909 |
|
361 |
plouden |
3941 |
\begin{equation}\label{tet_fit} |
362 |
plouden |
3942 |
q_{z} \approx q_{liq}+\frac{q_{ice}-q_{liq}}{2}\Bigg[\tanh\bigg(\frac{z-I_{L,m}}{w}\bigg)-\tanh\bigg(\frac{z-I_{R,m}}{w}\bigg)\Bigg]+\beta|(z-z_{mid})| |
363 |
gezelter |
3897 |
\end{equation} |
364 |
|
|
|
365 |
plouden |
3942 |
where $q_{liq}$ and $q_{ice}$ are fitting parameters for the local tetrahedral order parameter for the liquid and ice, $w$ is the width of the interface, and $I_{L,m}$ and $I_{R,m}$ are the midpoints of the left and right interface. The last term in \eqref{tet_fit} accounts for the influence the thermal gradient has on the tetrahedrality profile in the liquid region; where $\beta$ is a fitting parameter and $z_{mid}$ is the midpoint of the $z$-dimension of the simulation box. |
366 |
gezelter |
3897 |
|
367 |
plouden |
3942 |
In Figures \ref{fig:bComic} and \ref{fig:pComic} we see the $z$-dimensional profiles for several components of the basal and prismatic systems. In panel (a) of the figures we see the tetrahedrality profile of the systems (black circles). In the liquid region of the system, the local tetrahedral order parameter is approximately 0.75 while in the solid region the parameter is approximately 0.94, indicating a more tetrahedral structure of the water molecules. The hyperbolic tangent function used to fit the tetrahedrality profiles is in red and the verticle dotted lines denote the midpoint of the interfaces. The weak thermal gradient applied to the systems in order to keep the interface at a stable temperature, 225$\pm$5K, can be seen in panel (b). Lastly, the velocity gradient across the systems can be seen in panel (c). From panel (c), we can see liquid phase water molecules 10 \AA\ to 15 \AA\ from the interfaces are being dragged along with the ice block, indicating that the shearing of ice water is in the stick boundary condition. |
368 |
plouden |
3898 |
|
369 |
gezelter |
3914 |
\begin{figure} |
370 |
|
|
\includegraphics[width=\linewidth]{bComicStrip} |
371 |
|
|
\caption{\label{fig:bComic} The basal system: (a) The local tetrahedral order parameter,$q$, (black circles) fit by a hyperbolic tangent (red line), (b) the thermal gradient imposed on the system to maintain a stable interfacial temperature, and (c) the velocity gradient imposed on the system. The verticle dotted lines indicate the midpoint of the interfaces.} |
372 |
|
|
\end{figure} |
373 |
plouden |
3904 |
|
374 |
gezelter |
3914 |
\begin{figure} |
375 |
|
|
\includegraphics[width=\linewidth]{pComicStrip} |
376 |
|
|
\caption{\label{fig:pComic} The prismatic system: (a) The local tetrahedral order parameter,$q$, (black circles) fit by a hyperbolic tangent (red line), (b) the thermal gradient imposed on the system to maintain a stable interfacial temperature, and (c) the velocity gradient imposed on the system.The verticle dotted lines indicate the midpoint of the interfaces.} |
377 |
|
|
\end{figure} |
378 |
|
|
|
379 |
plouden |
3942 |
From the tetrahedrality fits, we found the interfacial width for the basal and prismatic systems to be 3.2$\pm$0.4 \AA\ and 3.6$\pm$0.2 \AA\ with no applied momentum flux. Over the range of shear rates investigated, 0.6$\pm$0.3 ms\textsuperscript{-1} to 5.3$\pm$0.5 ms\textsuperscript{-1} for the basal system and 0.9$\pm$0.2 ms\textsuperscript{-1} to 4.5$\pm$0.1 ms\textsuperscript{-1} for the prismatic, there was no appreciable change in the interface width found. The calculated values for the interfacial width over all shear rates investigated contained the control values within their error bars. |
380 |
gezelter |
3914 |
|
381 |
plouden |
3941 |
\subsubsection{Orientational Time Correlation Function} |
382 |
|
|
The orientational time correlation function (OTCF) gives insight of the local environment of molecules. The rate at which the function decays corresponds to how hindered the motions of a molecule are. The more hindered a molecules motion is the slower the function will decay, and the function decays more rapidly for molecules with less constrained motions. |
383 |
|
|
\begin{equation}\label{C(t)1} |
384 |
|
|
C_{2}(t)=\langle P_{2}(\mathbf{v}_{i}(t)\mathbf{v}_{i}(t=0))\rangle |
385 |
|
|
\end{equation} |
386 |
plouden |
3942 |
In eq. \eqref{C(t)1}, $P_{2}$ is the Legendre polynomial of the second order and $\mathbf{v}_{i}$ is the bisecting unit vector of the $i$th water molecule in the lab frame. |
387 |
gezelter |
3897 |
|
388 |
plouden |
3942 |
Here, we are evaluating this function across the $z$-dimension of the system as another measure of the change in the local environment and behavior of water molecules from the liquid region to the slushy interfacial region. After each of the 0.5 ns simulations, the systems were run for an additional 200 ps where the positions of every molecule in the system were recorded every 0.1 ps. The systems were then divided into 30 bins and the OTCF was evaluated for each bin. |
389 |
plouden |
3919 |
|
390 |
plouden |
3942 |
It has been shown that the OTCF for water can be fit by a triexponential decay\cite{Furse08}, where the three components of the decay correspond to a fast (<200 fs) reorientational piece driven by the restoring forces of existing hydrogen bonds, a middle (on the order of several ps) piece describing the large angle jumps that occur during the breaking and formation of new hydrogen bonds\cite{Laage08,Laage11}, and a slow (on the order of tens of ps) contribution describing the translational motion of the molecules. The OTCF data for each bin were truncated at 100 ps, and fit to the triexponential decay |
391 |
gezelter |
3897 |
\begin{equation} |
392 |
plouden |
3942 |
C_{2}(t) \approx a_{1}e^{-t/\tau_{short}}+a_{2}e^{-t/\tau_{middle}}+a_{3}e^{-t/\tau_{long}}+a_{4} |
393 |
gezelter |
3897 |
\end{equation} |
394 |
plouden |
3941 |
where $a_{1}+a_{2}+a_{3}+a_{4}=1$. An average value and standard deviation for each $\tau$ was obtained for each bin from the four runs. Lastly, the means and standard deviations were averaged about the center of the system. |
395 |
gezelter |
3897 |
|
396 |
plouden |
3941 |
\begin{figure} |
397 |
|
|
\includegraphics[width=\linewidth]{basal_Tau_comic_strip} |
398 |
|
|
\caption{\label{fig:basal_Tau_comic_strip} The orientational time correlation function for the basal system fit by the triexponential decay $a_{1}e^{-t/\tau_{short}}+a_{2}e^{-t/\tau_{middle}}+a_{3}e^{-t/\tau_{long}}+a_{4}$ where $a_{1}+a_{2}+a_{3}+a_{4}=1$. The verticle dotted line indicates the average midpoint of the interface as determined by the tetrahedrality fit. In (b) and (c) $\tau_{middle}$ and $\tau_{long}$ have a consistent value of about 5.5 ps and 50 ps in the liquid region, and increase in value approaching the interface. $\tau_{middle}$ corresponds to the breaking and making of hydrogen bonds as explained by extended jump model proposed by Laage and Hynes\cite{Laage08,Laage11} and $\tau_{long}$ relates to the translational motion of the molecules. In (a), we see that $\tau_{short}$ has a value of about 71 fs in the liquid region, and decreases in value approaching the interface. This component of the decay corresponds to the reorientational forces of existing hydrogen bonds on the inertial rotational of the molecules. Thus $\tau_{short}$ decreases approaching the interface due to the hindered range of motion of the molecules. } |
399 |
|
|
\end{figure} |
400 |
gezelter |
3897 |
|
401 |
gezelter |
3914 |
\begin{figure} |
402 |
plouden |
3941 |
\includegraphics[width=\linewidth]{prismatic_Tau_comic_strip} |
403 |
|
|
\caption{\label{fig:prismatic_Tau_comic_strip} The orientational time correlation function for the prismatic system fit by the triexponential decay $a_{1}e^{-t/\tau_{short}}+a_{2}e^{-t/\tau_{middle}}+a_{3}e^{-t/\tau_{long}}+a_{4}$ where $a_{1}+a_{2}+a_{3}+a_{4}=1$. The verticle dotted line indicates the average midpoint of the interface as determined by the tetrahedrality fit. In (b) and (c) $\tau_{middle}$ and $\tau_{long}$ have a consistent value of about 3.5 ps and 30 ps in the liquid region, and increase in value approaching the interface. $\tau_{middle}$ corresponds to the breaking and making of hydrogen bonds as explained by extended jump model proposed by Laage and Hynes\cite{Laage08,Laage11} and $\tau_{long}$ relates to the translational motion of the molecules. In (a), we see that $\tau_{short}$ has a value of about 73 fs in the liquid region, and decreases in value approaching the interface. This component of the decay corresponds to the reorientational forces of existing hydrogen bonds on the inertial rotational of the molecules. Thus $\tau_{short}$ decreases approaching the interface due to the hindered range of motion of the molecules.} |
404 |
gezelter |
3914 |
\end{figure} |
405 |
plouden |
3904 |
|
406 |
plouden |
3942 |
Figures \ref{fig:basal_Tau_comic_strip} and \ref{fig:prismatic_Tau_comic_strip} plot the decomposition of the OTCF at varying displacements from the center of the ice for the basal and prismatic systems. We see in (a) $\tau_{short}$, (b) $\tau_{middle}$, and (c) $\tau_{long}$ for the control system (no applied momentum flux) in black, and a system with a large shear rate in red. The verticle dotted lines at a displacement of about 17 \AA\ and 9 \AA\ denote the midpoints of the interfaces as determined by the hyperbolic tangent fit of the tetrahedrality profile. |
407 |
plouden |
3937 |
|
408 |
plouden |
3942 |
In panels (a), we see at large displacements from the center of the ice $\tau_{short}$ for the basal system has a value of about 71 fs and 72 fs for the prismatic. Decreasing in displacement from about 26 \AA\ to about 19 \AA\ in the basal system, the value of $\tau_{short}$ decreases to about 63 fs. Likewise, $\tau_{short}$ decreases to about 63 fs from roughly 20 \AA\ to 12 \AA\ . This is due to the increasingly constrained motion of the water molecules as we approach the interface. In panels (b), $\tau_{middle}$ at large displacements from the ice has a value of about 5.5 ps and 3 ps for the basal and prismatic systems. We find $\tau_{middle}$ increases in value as we approach the interface in both cases. This component of the decay corresponds to the rearrangement of the hydrogen bonding network, which takes longer as the molecules motion becomes more constrained. In panels (c), $\tau_{long}$ has a value of about 50 ps for the basal system and roughly 30 ps for the prismatic system at large displacements from the interface. Similar to $\tau_{middle}$, $\tau_{long}$ also increases in value as we approach the interface for both systems. It is also apparent from Figures \ref{fig:basal_Tau_comic_strip} and \ref{fig:prismatic_Tau_comic_strip} that shearing the ice water has no effect on the orientational decay time, or on any of its decomposed components. |
409 |
plouden |
3937 |
|
410 |
plouden |
3942 |
For each system, there is an apparent approximate value for $\tau_{short}$, $\tau_{middle}$, and $\tau_{long}$ at large displacements from the interface. There also appears to be a single displacement, $d_{basal}$ and $d_{prismatic}$, from the interface at which all three decay times begin to deviate from their bulk liquid values. We found $d_{basal}$ and $d_{prismatic}$ to be roughly 15 \AA\ and 8 \AA\ respectively. These two results indicate that the dynamics of the water molecules within $d_{basal}$ and $d_{prismatic}$ are being significantly perturbed by the ice and/or the interface, even though the structural width of the interface by analysis of the tetrahedrality profile indicates that bulk liquid structure of water is recovered after about 4 \AA\ from the edge of the ice. |
411 |
plouden |
3937 |
|
412 |
plouden |
3942 |
Beaglehole and Wilson have measured the ice/water interface to have a thickness approximately 10 \AA\ for both the basal and prismatic face of ice by ellipticity measurements \cite{Beaglehole93}. Structurally, we have found the basal and prismatic interfacial width to be 3.2$\pm$0.4 \AA\ and 3.6$\pm$0.2 \AA\ . However, we have shown through decomposition of the OTCF a much larger interfacial region exists in which the dynamics of the water molecules behave differently than those of the bulk liquid. |
413 |
plouden |
3937 |
|
414 |
plouden |
3941 |
\subsection{Coefficient of Friction of the Interface} |
415 |
plouden |
3942 |
As the ice is sheared through the liquid, there will be a friction between the solid and the liquid. Pit has shown how to calculate the coefficient of friction $\lambda$ for a solid-liquid interface for a Newtonian fluid of viscosity $\eta$ and has a slip length of $\delta$. \cite{Pit99} |
416 |
|
|
\begin{equation}\label{Pit} |
417 |
|
|
\lambda=\eta/\delta |
418 |
|
|
\end{equation} |
419 |
|
|
From linear response theory, $\eta$ can be obtained from the imposed momentum flux and the slope of the velocity about the dimension of the imposed flux.\cite{Kuang12} |
420 |
|
|
\begin{equation}\label{Kuang} |
421 |
|
|
j_{z}(p_{x})=-\eta\frac{\partial v_{x}}{\partial z} |
422 |
|
|
\end{equation} |
423 |
|
|
Solving eq. \eqref{Kuang} for $\eta$ and substituting the result into eq. \eqref{Pit}, we obtain an alternate expression for the coefficient of friction. |
424 |
plouden |
3941 |
\begin{equation} |
425 |
plouden |
3942 |
\lambda=-\frac{j_{z}(p_{x})}{\delta \frac{\partial v_{x}}{\partial z}} |
426 |
plouden |
3941 |
\end{equation} |
427 |
plouden |
3937 |
|
428 |
plouden |
3942 |
For our simulations, we obtain $\delta$ from the difference between the structural edge of the ice block determined by the tetrahedrality profile fit, and the intersection of the linear regression of the $v_{x}$ profiles about the $z$-dimension for the ice and liquid. (See Figure \ref{fig:delta_example}) The coefficient of friction for the basal and the prismatic facets were found to be (0.07$\pm$0.01) amu fs\textsuperscript{-1} and (0.032$\pm$0.007) amu fs\textsuperscript{-1}. It is known that the basal and prismatic faces have different surface structures. The basal face is smoother than the prismatic with small alternating valleys and crests, while the prismatic surface has deep corrugating channels. We believe the reason that the prismatic face's coefficient of friction was found to be smaller than the basal's is due to the direction of the shear. The shear of the ice/water was in the same direction of the corrugating channels, allowing water molecules to pass through the channels during the shear. |
429 |
plouden |
3939 |
|
430 |
plouden |
3942 |
\begin{figure} |
431 |
|
|
\includegraphics[width=\linewidth]{delta_example} |
432 |
|
|
\caption{\label{fig:delta_example} A schematic of determining the slip length ($\delta$). The slip length is the difference of the structural starting point of the ice and the point of intersection of the linear regressions of the liquid phase velocity profile (red) and of the solid ice velocity profile (black). The dotted line indicates the location of the ice as determined by the tetrahedrality profile.} |
433 |
|
|
\end{figure} |
434 |
|
|
|
435 |
|
|
|
436 |
gezelter |
3897 |
\section{Conclusion} |
437 |
gezelter |
3945 |
Here we have simulated the basal and prismatic facets of an SPC/E model of the ice I$_\mathrm{h}$ / water interface. Using VSS-RNEMD, the ice was sheared relative to the liquid while imposed thermal gradients kept the interface at a stable temperature. Caculation of the local tetrahedrality order parameter has shown an apparent independence of the shear rate on the interfacial width, which was found to be 3.2$\pm$0.4 \AA\ and 3.6$\pm$0.2 \AA\ for the basal and prismatic systems. The orientational time correlation function was calculated from varying displacements from the interface. Decomposition by a triexponential decay also showed an apparent independence of the shear rate. The short time decay due to the restoring forces of existing hydrogen bonds decreased at close displacements from the interface, while the middle and long time decays were found to increase. There is also an apparent displacement, $d_{basal}$ and $d_{prismatic}$, from the interface at which these deviations from bulk liquid values occurs. We found $d_{basal}$ and $d_{prismatic}$ to be approximately 15 \AA\ and 8 \AA\ . This implies that the dynamics of water molecules which are structurally equivalent to bulk phase molecules are being perturbed by the presence of the ice and/or the interface. The coefficient of friction of each of the facets was also determined. They were found to be (0.07$\pm$0.01) amu fs\textsuperscript{-1} and (0.032$\pm$0.007) amu fs\textsuperscript{-1} for the basal and prismatic facets respectively. We believe the large difference between the two friction coefficients is due to the direction of the shear and the surface structure of the crystal facets. |
438 |
gezelter |
3897 |
|
439 |
gezelter |
3914 |
\begin{acknowledgement} |
440 |
|
|
Support for this project was provided by the National Science |
441 |
|
|
Foundation under grant CHE-0848243. Computational time was provided |
442 |
|
|
by the Center for Research Computing (CRC) at the University of |
443 |
|
|
Notre Dame. |
444 |
|
|
\end{acknowledgement} |
445 |
gezelter |
3897 |
|
446 |
gezelter |
3914 |
\newpage |
447 |
|
|
\bibstyle{achemso} |
448 |
plouden |
3909 |
\bibliography{iceWater} |
449 |
gezelter |
3897 |
|
450 |
gezelter |
3914 |
\begin{tocentry} |
451 |
|
|
\begin{wrapfigure}{l}{0.5\textwidth} |
452 |
|
|
\begin{center} |
453 |
|
|
\includegraphics[width=\linewidth]{SystemImage.png} |
454 |
|
|
\end{center} |
455 |
|
|
\end{wrapfigure} |
456 |
|
|
An image of our system. |
457 |
|
|
\end{tocentry} |
458 |
gezelter |
3897 |
|
459 |
gezelter |
3914 |
\end{document} |
460 |
gezelter |
3897 |
|
461 |
plouden |
3924 |
%************************************************************** |
462 |
|
|
%Non-mass weighted slopes (amu*Angstroms^-2 * fs^-1) |
463 |
|
|
% basal: slope=0.090677616, error in slope = 0.003691743 |
464 |
|
|
%prismatic: slope = 0.050101506, error in slope = 0.001348181 |
465 |
|
|
%Mass weighted slopes (Angstroms^-2 * fs^-1) |
466 |
|
|
%basal slope = 4.76598E-06, error in slope = 1.94037E-07 |
467 |
|
|
%prismatic slope = 3.23131E-06, error in slope = 8.69514E-08 |
468 |
|
|
%************************************************************** |