trunk/iceWater2/iceWater2.tex

% ****** Start of file aipsamp.tex ******
%
%   This file is part of the AIP files in the AIP distribution for REVTeX 4.
%   Version 4.1 of REVTeX, October 2009
%
%   Copyright (c) 2009 American Institute of Physics.
%
%   See the AIP README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.1
%
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex  aipsamp
%  2)  bibtex aipsamp
%  3)  latex  aipsamp
%  4)  latex  aipsamp
%
% Use this file as a source of example code for your aip document.
% Use the file aiptemplate.tex as a template for your document.
\documentclass[%
 aip,jcp,
 amsmath,amssymb,
preprint,%
% reprint,%
%author-year,%
%author-numerical,%
jcp]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
%\usepackage{bm}% bold math
\usepackage{times}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}

\begin{document}

\title{Friction at Water / Ice-I$_\mathrm{h}$ interfaces: Do the
  Facets of Ice Have Different Hydrophilicity?}

\author{Patrick B. Louden}

\author{J. Daniel Gezelter}
 \email{gezelter@nd.edu.}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}

\date{\today}

\begin{abstract}
Abstract abstract abstract...
\end{abstract}

\maketitle

\section{Introduction}
Explain a little bit about ice Ih, point group stuff. 

Mention previous work done / on going work by other people. Haymet and Rick
seem to be investigating how the interfaces is perturbed by the presence of 
ions. This is the conlcusion of a recent publication of the basal and 
prismatic facets of ice Ih, now presenting the pyramidal and secondary 
prism facets under shear.

\section{Methodology}

\begin{figure}
\includegraphics[width=\linewidth]{SP_comic_strip}
\caption{\label{fig:spComic} The secondary prism interface with a shear
 rate of 3.5 ms\textsuperscript{-1}. Lower panel: the local tetrahedral order
parameter, $q(z)$, (black circles) and the hyperbolic tangent fit (red line).
Middle panel: the imposed thermal gradient required to maintain a fixed
interfacial temperature. Upper panel: the transverse velocity gradient that
develops in response to an imposed momentum flux. The vertical dotted lines
indicate the locations of the midpoints of the two interfaces.}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{Pyr_comic_strip}
\caption{\label{fig:pyrComic} The pyramidal interface with a shear rate of 3.8 \
ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:spComic}.}
\end{figure}

\subsection{Pyramidal and secondary prism system construction}

The construction of the pyramidal and secondary prism systems follows that of 
the basal and prismatic systems presented elsewhere\cite{Louden13}, however
the ice crystals and water boxes were equilibrated and combined at 50K 
instead of 225K. The ice / water systems generated were then equilibrated 
to 225K. The resulting pyramidal system was 
$37.47 \times 29.50 \times 93.02$ \AA\ with 1216
SPC/E molecules in the ice slab, and 2203 in the liquid phase. The secondary 
prism system generated was $71.87 \times 31.66 \times 161.55$ \AA\ with 3840 
SPC/E molecules in the ice slab and 8176 molecules in the liquid phase.

\subsection{Computational details}
% Do we need to justify the sims at 225K? 
% No crystal growth or shrinkage over 2 successive 1 ns NVT simulations for
%    either the pyramidal or sec. prism ice/water systems.

The computational details performed here were equivalent to those reported
in the previous publication\cite{Louden13}. The only changes made to the 
previously reported procedure were the following. VSS-RNEMD moves were 
attempted every 2 fs instead of every 50 fs. This was done to minimize
the magnitude of each individual VSS-RNEMD perturbation to the system.

All pyramidal simulations were performed under the NVT ensamble except those
during which statistics were accumulated for the orientational correlation
function, which were performed under the NVE ensamble. All secondary prism 
simulations were performed under the NVE ensamble. 

\section{Results and discussion}
\subsection{Interfacial width}
In the literature there is good agreement that between the solid ice and 
the bulk water, there exists a region of 'slush-like' water molecules. 
In this region, the water molecules are structurely distinguishable and 
behave differently than those of the solid ice or the bulk water.
The characteristics of this region have been defined by both structural
and dynamic properties; and its width has been measured by the change of these 
properties from their bulk liquid values to those of the solid ice. 
Examples of these properties include the density, the diffusion constant, and 
the translational order profile. \cite{Bryk02,Karim90,Gay02,Hayward01,Hayward02,Karim88}  

Since the VSS-RNEMD moves used to impose the thermal and velocity gradients
 perturb the momenta of the water molecules in 
the systems, parameters that depend on translational motion may give
faulty results. A stuructural parameter will be less effected by the 
VSS-RNEMD perturbations to the system. Due to this, we have used the 
local order tetrahedral parameter to quantify the width of the interface,
 which was originally described by Kumar\cite{Kumar09} and 
Errington\cite{Errington01} and explained in our
previous publication\cite{Louden13} in relation to an ice/water system.

Paragraph and eq. for tetrahedrality here. 

To determine the width of the interfaces, each of the systems were 
divided into 100 artificial bins along the 
$z$-dimension, and the local tetrahedral order parameter, $q(z)$, was 
time-averaged for each of the bins, resulting in a tetrahedrality profile of 
the system. These profiles are shown across the $z$-dimension of the systems
in panel $a$ of Figures \ref{fig:spComic}
and \ref{fig:pyrComic} (black circles). The $q(z)$ function has a range of 
(0,1), where a larger value indicates a more tetrahedral environment.
The $q(z)$ for the bulk liquid was found to be $\approx $ 0.77, while values of
$\approx $0.92 were more common for the ice. The tetrahedrality profiles were
fit using a hyperbolic tangent\cite{Louden13} designed to smoothly fit the 
bulk to ice 
transition, while accounting for the thermal influence on the profile by the
kinetic energy exchanges of the VSS-RNEMD moves. In panels $b$ and $c$, the
imposed thermal and velocity gradients can be seen. The verticle dotted
lines traversing all three panels indicate the midpoints of the interface
as determined by the hyperbolic tangent fit of the tetrahedrality profiles.
 
From fitting the tetrahedrality profiles for each of the 0.5 nanosecond 
simulations (panel c of Figures \ref{fig:spComic} and \ref{fig:pyrComic}) 
by Eq. 6\cite{Louden13},we find the interfacial width to be
 $3.2 \pm 0.2$ and $3.2 \pm 0.2$ \AA\ for the control system with no applied 
momentum flux for both the pyramidal and secondary prism systems. 
Over the range of shear rates investigated, 
$0.6 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 5.6 \pm 0.4 \mathrm{ms}^{-1}$ for 
the pyramidal system and $0.9 \pm 0.3 \mathrm{ms}^{-1} \rightarrow 5.4 \pm 0.1
\mathrm{ms}^{-1}$ for the secondary prism, we found no significant change in
the interfacial width. This follows our previous findings of the basal and
prismatic systems, in which the interfacial width was invarient of the
shear rate of the ice. The interfacial width of the quiescent basal and 
prismatic systems was found to be $3.2 \pm 0.4$ \AA\ and $3.6 \pm 0.2$ \AA\ 
respectively. Over the range of shear rates investigated, $0.6 \pm 0.3 
\mathrm{ms}^{-1} \rightarrow 5.3 \pm 0.5 \mathrm{ms}^{-1}$ for the basal 
system and $0.9 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 4.5 \pm 0.1 
\mathrm{ms}^{-1}$ for the prismatic.

These results indicate that the surface structure of the exposed ice crystal
has little to no effect on how far into the bulk the ice-like structural 
ordering is. Also, it appears that the interface is not structurally effected
by shearing the ice through water. 


\subsection{Orientational dynamics}
%Should we include the math here?
The orientational time correlation function,
\begin{equation}\label{C(t)1}
  C_{2}(t)=\langle P_{2}(\mathbf{u}(0)\cdot \mathbf{u}(t))\rangle,
\end{equation}
helps indicate the local environment around the water molecules. The function 
begins with an initial value of unity, and decays to zero as the water molecule
loses memory of its former orientation. Observing the rate at which this decay
occurs can provide insight to the mechanism and timescales for the relaxation.
In eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial, and 
$\mathbf{u}$ is the the bisecting HOH vector. The angle brackets indicate
an ensemble average over all the water molecules in a given spatial region.

To investigate the dynamics of the water molecules across the interface, the 
systems were divided in the $z$-dimension into bins, each $\approx$ 3 \AA\
 wide, and \eqref{C(t)1} was computed for each of the bins. A water 
molecule was allocated to a particular bin if it was initially in the bin
at time zero. To compute \eqref{C(t)1}, each 0.5 ns simulation was followed 
by an additional 200 ps microcanonical (NVE) simulation during which the 
position and orientations of each molecule were recorded every 0.1 ps.
 
The data obtained for each bin was then fit to a triexponential decay given by
\begin{equation}\label{C(t)_fit}
C_{2}(t) \approx a e^{-t/\tau_\mathrm{short}} + b e^{-t/\tau_\mathrm{middle}} +\
 c
e^{-t/\tau_\mathrm{long}} + (1-a-b-c)
\end{equation}
where $\tau_{short}$ corresponds to the librational motion of the water  
molecules, $\tau_{middle}$ corresponds to jumps between the breaking and 
making of hydrogen bonds, and $\tau_{long}$ corresponds to the translational
motion of the water molecules. The last term in \eqref{C(t)_fit} accounts
for the water molecules trapped in the ice which do not experience any
long-time orientational decay.

In Figures \ref{fig:PyrOrient} and \ref{fig:SPorient} we see the $z$-coordinate
profiles for the three decay constants, $\tau_{short}$ (panel a), 
$\tau_{middle}$ (panel b),
and $\tau_{long}$ (panel c) for the pyramidal and secondary prismatic systems
respectively. The control experiments (no shear) are shown in black, and 
an experiment with an imposed momentum flux is shown in red. The vertical
dotted line traversing all three panels denotes the midpoint of the 
interface as determined by the local tetrahedral order parameter fitting.
In the liquid regions of both systems, we see that $\tau_{middle}$ and
$\tau_{long}$ have approximately consistent values of $3-6$ ps and $30-40$ ps,
resepctively, and increase in value as we approach the interface. Conversely,
in panel a, we see that $\tau_{short}$ decreases from the liquid value 
of $72-76$ fs as we approach the interface. We believe this speed up is due to
the constrained motion of librations closer to the interface. Both the 
approximate values for the decays and relative trends  match those reported 
previously for the basal and prismatic interfaces.

As done previously, we have attempted to quantify the distance, $d_{pyramidal}$
and $d_{secondary prism}$, from the 
interface that the deviations from the bulk liquid values begin. This was done
by fitting the orientational decay constant $z$-profiles by
\begin{equation}\label{tauFit}
\tau(z)\approx\tau_{liquid}+(\tau_{solid}-\tau_{liquid})e^{-(z-z_{wall})/d}
\end{equation}
where $\tau_{liquid}$ and $\tau_{solid}$ are the liquid and projected solid
values of the decay constants, $z_{wall}$ is the location of the interface,
and $d$ is the displacement from the interface at which these deviations
occur. The values for $d_{pyramidal}$ and $d_{secondary prismatic}$ were 
determined
for each of the decay constants, and then averaged for better statistics 
($\tau_{middle}$ was ommitted for secondary prism). For the pyramidal system,
$d_{pyramidal}$ was found to be 2.7 \AA\ for both the control and the sheared
system. We found $d_{secondary prismatic}$ to be slightly larger than 
$d_{pyramidal}$ for both the control and with an applied shear, with 
displacements of $4$ \AA\ for the control system and $3$ \AA\ for the 
experiment with the imposed momentum flux. These values are consistent with
those found for the basal ($d_{basal}\approx2.9$ \AA\ ) and prismatic 
($d_{prismatic}\approx3.5$ \AA\ ) systems.

\subsection{Coefficient of friction of the interfaces}
While investigating the kinetic coefficient of friction for the larger
prismatic system, there was found to be a dependence for $\mu_k$ 
on the temperature of the liquid water in the system. We believe this 
dependence
arrises from the sharp discontinuity of the viscosity for the SPC/E model
at temperatures approaching 200 K\cite{kuang12}. Due to this, we propose
a weighting to the structural interfacial parameter, $\kappa$ by the 
viscosity at $225$ K, the temperature of the interface. $\kappa$ is 
traditionally defined as
\begin{equation}\label{kappa}
\kappa = \eta/\delta
\end{equation}
where $\eta$ is the viscosity and $\delta$ is the slip length. 
In our ice/water shearing simulations, the system has reached a steady state 
when the applied force,

\begin{equation}
f_{applied} = \mathbf{j}_z(\mathbf{p})L_x L_y
\end{equation} 
is equal to the frictional force resisting the motion of the ice block
\begin{equation}
f_{friction} = \frac{\eta \mathbf{v} L_x L_y}{\delta} 
\end{equation}
where $\mathbf{v}$ is the relative velocity of the liquid from the ice.
When this condition is met, we are able to solve the resulting expression to 
obtain,
\begin{equation}\label{force_equality}
\frac{\eta}{\delta} = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}}
\end{equation}
From \eqref{kappa}, \eqref{force_equality} becomes
\begin{equation}
\kappa = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}}
\end{equation}
which we will multiply by a viscosity weighting term to reach
\begin{equation} \label{kappa2}
\kappa = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}} \frac{\eta(225)}{\eta(T)}
\end{equation}
Assuming linear response theory is valid, an expression for ($\eta$) can
be found from the imposed momentum flux and the measured velocity gradient.
\begin{equation}\label{eta_eq}
\eta = \frac{\mathbf{j}_z(\mathbf{p})}{\frac{\partial v_x}{\partial z}}
\end{equation}
Substituting eq \eqref{eta_eq} into eq \eqref{kappa2} we arrive at
\begin{equation}
  \kappa = \frac{\frac{\partial v_x}{\partial z}}{\mathbf{v}}\eta(225)
\end{equation} 

\begin{table}[h]
\centering
\caption{$\kappa$ values for the basal, prismatic, pyramidal, and secondary    prismatic facets of Ice-I$_\mathrm{h}$}
\label{tab:kappa}
\begin{tabular}{|ccc|}  \hline
           & \multicolumn{2}{c|}{$\kappa_{Drag direction}$} \\
 Interface & $\kappa_{x}$     & $\kappa_{y}$   \\ \hline
     basal & $0.00059 \pm 0.00003$ & $0.00065 \pm 0.00008$ \\
 prismatic & $0.00030 \pm 0.00002$ & $0.00030 \pm 0.00001$ \\
 pyramidal & $0.00058 \pm 0.00004$ & $0.00061 \pm 0.00005$ \\
 secondary prism & $0.00035 \pm 0.00001$ & $0.00033 \pm 0.00002$ \\ \hline
\end{tabular}
\end{table}


To obtain the value of $\eta(225)$ for the SPC/E model, a $31.09 \times 29.38
\times 124.39$ \AA\ box with 3744 water molecules was equilibrated to 225K, 
and 5 unique shearing experiments were performed. Each experiment was
conducted in the microcanonical ensemble (NVE) and were 5 ns in 
length. The VSS were attempted every timestep, which was set to 2 fs. 
For our SPC/E systems, we found $\eta(225)$  to be $0.0148 \pm 0.0007 Pa s$,
roughly ten times larger than the value found for 280 K SPC/E water by
Kuang\cite{kuang12}. 


\begin{table}[h]
\centering
\caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript\
{-1}). \\ 
\textsuperscript{a} See ref. \onlinecite{Louden13}. }
\label{tab:lambda}
\begin{tabular}{|ccc|}  \hline
           & \multicolumn{2}{c|}{Drag direction} \\
 Interface & $x$               & $y$  \\ \hline
     basal\textsuperscript{a} &  $0.08 \pm 0.02$  & $0.09 \pm 0.03$ \\
 prismatic (T = 225)\textsuperscript{a} & $0.037 \pm 0.008$ & $0.04 \pm 0.01$ \\
 prismatic (T = 230) & $0.10 \pm 0.01$ & $0.070 \pm 0.006$\\
 pyramidal & $0.13 \pm 0.03$   & $0.14 \pm 0.03$ \\
 secondary prism & $0.13 \pm 0.02$ & $0.12 \pm 0.03$ \\ \hline
\end{tabular}
\end{table}
 

\begin{figure}
\includegraphics[width=\linewidth]{Pyr-orient}
\caption{\label{fig:PyrOrient} The three decay constants of the 
orientational time correlation function, $C_2(t)$, for water as a function
of distance from the center of the ice slab. The vertical dashed line
indicates the edge of the pyramidal ice slab determined by the local order
tetrahedral parameter. The control (black circles) and sheared (red squares) 
experiments were fit by a shifted exponential decay (Eq. 9\cite{Louden13})
shown by the black and red lines respectively. The upper two panels show that 
translational and hydrogen bond making and breaking events slow down 
through the interface while approaching the ice slab. The bottom most panel
shows the librational motion of the water molecules speeding up approaching
the ice block due to the confined region of space allowed for the molecules
 to move in.}
\end{figure}  

\begin{figure}
\includegraphics[width=\linewidth]{SP-orient-less}
\caption{\label{fig:SPorient} Decay constants for $C_2(t)$ at the secondary 
prism face. Panel descriptions match those in \ref{fig:PyrOrient}.}
\end{figure}


\section{Conclusion}
Conclude conclude conclude...


\begin{acknowledgments}
  Support for this project was provided by the National
  Science Foundation under grant CHE-1362211. Computational time was
  provided by the Center for Research Computing (CRC) at the
  University of Notre Dame.
\end{acknowledgments}

\newpage

\bibliography{iceWater}

\end{document}
Revision:	4219
Committed:	Fri Sep 12 20:59:54 2014 UTC (10 years, 10 months ago) by plouden
Content type:	application/x-tex
File size:	18074 byte(s)
Log Message:	added kappa table and other changes to iceWater2.tex
#	User	Rev	Content
1	gezelter	4217	% **** Start of file aipsamp.tex ****
2			%
3			% This file is part of the AIP files in the AIP distribution for REVTeX 4.
4			% Version 4.1 of REVTeX, October 2009
5			%
6			% Copyright (c) 2009 American Institute of Physics.
7			%
8			% See the AIP README file for restrictions and more information.
9			%
10			% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
11			% as well as the rest of the prerequisites for REVTeX 4.1
12			%
13			% It also requires running BibTeX. The commands are as follows:
14			%
15			% 1) latex aipsamp
16			% 2) bibtex aipsamp
17			% 3) latex aipsamp
18			% 4) latex aipsamp
19			%
20			% Use this file as a source of example code for your aip document.
21			% Use the file aiptemplate.tex as a template for your document.
22			\documentclass[%
23			aip,jcp,
24			amsmath,amssymb,
25			preprint,%
26			% reprint,%
27			%author-year,%
28			%author-numerical,%
29			jcp]{revtex4-1}
30	plouden	4192
31	gezelter	4217	\usepackage{graphicx}% Include figure files
32			\usepackage{dcolumn}% Align table columns on decimal point
33			%\usepackage{bm}% bold math
34			\usepackage{times}
35			\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions
36	plouden	4192	\usepackage{url}
37
38			\begin{document}
39
40	gezelter	4217	\title{Friction at Water / Ice-I$_\mathrm{h}$ interfaces: Do the
41			Facets of Ice Have Different Hydrophilicity?}
42	plouden	4192
43	gezelter	4217	\author{Patrick B. Louden}
44	plouden	4192
45	gezelter	4217	\author{J. Daniel Gezelter}
46			\email{gezelter@nd.edu.}
47			\affiliation{Department of Chemistry and Biochemistry, University
48			of Notre Dame, Notre Dame, IN 46556}
49
50	plouden	4192	\date{\today}
51
52			\begin{abstract}
53			Abstract abstract abstract...
54			\end{abstract}
55
56	gezelter	4217	\maketitle
57	plouden	4192
58			\section{Introduction}
59			Explain a little bit about ice Ih, point group stuff.
60
61			Mention previous work done / on going work by other people. Haymet and Rick
62			seem to be investigating how the interfaces is perturbed by the presence of
63			ions. This is the conlcusion of a recent publication of the basal and
64			prismatic facets of ice Ih, now presenting the pyramidal and secondary
65			prism facets under shear.
66
67			\section{Methodology}
68
69			\begin{figure}
70			\includegraphics[width=\linewidth]{SP_comic_strip}
71			\caption{\label{fig:spComic} The secondary prism interface with a shear
72			rate of 3.5 ms\textsuperscript{-1}. Lower panel: the local tetrahedral order
73			parameter, $q(z)$, (black circles) and the hyperbolic tangent fit (red line).
74			Middle panel: the imposed thermal gradient required to maintain a fixed
75			interfacial temperature. Upper panel: the transverse velocity gradient that
76			develops in response to an imposed momentum flux. The vertical dotted lines
77			indicate the locations of the midpoints of the two interfaces.}
78			\end{figure}
79
80			\begin{figure}
81			\includegraphics[width=\linewidth]{Pyr_comic_strip}
82			\caption{\label{fig:pyrComic} The pyramidal interface with a shear rate of 3.8 \
83			ms\textsuperscript{-1}. Panel descriptions match those in figure \ref{fig:spComic}.}
84			\end{figure}
85
86			\subsection{Pyramidal and secondary prism system construction}
87
88			The construction of the pyramidal and secondary prism systems follows that of
89			the basal and prismatic systems presented elsewhere\cite{Louden13}, however
90	plouden	4194	the ice crystals and water boxes were equilibrated and combined at 50K
91			instead of 225K. The ice / water systems generated were then equilibrated
92			to 225K. The resulting pyramidal system was
93	plouden	4192	$37.47 \times 29.50 \times 93.02$ \AA\ with 1216
94			SPC/E molecules in the ice slab, and 2203 in the liquid phase. The secondary
95			prism system generated was $71.87 \times 31.66 \times 161.55$ \AA\ with 3840
96			SPC/E molecules in the ice slab and 8176 molecules in the liquid phase.
97
98			\subsection{Computational details}
99			% Do we need to justify the sims at 225K?
100			% No crystal growth or shrinkage over 2 successive 1 ns NVT simulations for
101			% either the pyramidal or sec. prism ice/water systems.
102
103			The computational details performed here were equivalent to those reported
104			in the previous publication\cite{Louden13}. The only changes made to the
105			previously reported procedure were the following. VSS-RNEMD moves were
106	plouden	4194	attempted every 2 fs instead of every 50 fs. This was done to minimize
107			the magnitude of each individual VSS-RNEMD perturbation to the system.
108	plouden	4192
109			All pyramidal simulations were performed under the NVT ensamble except those
110			during which statistics were accumulated for the orientational correlation
111			function, which were performed under the NVE ensamble. All secondary prism
112			simulations were performed under the NVE ensamble.
113
114			\section{Results and discussion}
115	plouden	4194	\subsection{Interfacial width}
116			In the literature there is good agreement that between the solid ice and
117			the bulk water, there exists a region of 'slush-like' water molecules.
118	plouden	4219	In this region, the water molecules are structurely distinguishable and
119	plouden	4194	behave differently than those of the solid ice or the bulk water.
120			The characteristics of this region have been defined by both structural
121	plouden	4215	and dynamic properties; and its width has been measured by the change of these
122	plouden	4194	properties from their bulk liquid values to those of the solid ice.
123			Examples of these properties include the density, the diffusion constant, and
124	gezelter	4217	the translational order profile. \cite{Bryk02,Karim90,Gay02,Hayward01,Hayward02,Karim88}
125	plouden	4192
126	plouden	4215	Since the VSS-RNEMD moves used to impose the thermal and velocity gradients
127			perturb the momenta of the water molecules in
128			the systems, parameters that depend on translational motion may give
129	plouden	4194	faulty results. A stuructural parameter will be less effected by the
130	plouden	4215	VSS-RNEMD perturbations to the system. Due to this, we have used the
131			local order tetrahedral parameter to quantify the width of the interface,
132			which was originally described by Kumar\cite{Kumar09} and
133			Errington\cite{Errington01} and explained in our
134			previous publication\cite{Louden13} in relation to an ice/water system.
135	plouden	4194
136	plouden	4215	Paragraph and eq. for tetrahedrality here.
137
138			To determine the width of the interfaces, each of the systems were
139			divided into 100 artificial bins along the
140	plouden	4194	$z$-dimension, and the local tetrahedral order parameter, $q(z)$, was
141			time-averaged for each of the bins, resulting in a tetrahedrality profile of
142			the system. These profiles are shown across the $z$-dimension of the systems
143			in panel $a$ of Figures \ref{fig:spComic}
144			and \ref{fig:pyrComic} (black circles). The $q(z)$ function has a range of
145			(0,1), where a larger value indicates a more tetrahedral environment.
146	plouden	4215	The $q(z)$ for the bulk liquid was found to be $\approx $ 0.77, while values of
147	plouden	4194	$\approx $0.92 were more common for the ice. The tetrahedrality profiles were
148			fit using a hyperbolic tangent\cite{Louden13} designed to smoothly fit the
149			bulk to ice
150			transition, while accounting for the thermal influence on the profile by the
151			kinetic energy exchanges of the VSS-RNEMD moves. In panels $b$ and $c$, the
152			imposed thermal and velocity gradients can be seen. The verticle dotted
153			lines traversing all three panels indicate the midpoints of the interface
154			as determined by the hyperbolic tangent fit of the tetrahedrality profiles.
155
156	plouden	4192	From fitting the tetrahedrality profiles for each of the 0.5 nanosecond
157	plouden	4194	simulations (panel c of Figures \ref{fig:spComic} and \ref{fig:pyrComic})
158	plouden	4215	by Eq. 6\cite{Louden13},we find the interfacial width to be
159			$3.2 \pm 0.2$ and $3.2 \pm 0.2$ \AA\ for the control system with no applied
160			momentum flux for both the pyramidal and secondary prism systems.
161			Over the range of shear rates investigated,
162	plouden	4192	$0.6 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 5.6 \pm 0.4 \mathrm{ms}^{-1}$ for
163			the pyramidal system and $0.9 \pm 0.3 \mathrm{ms}^{-1} \rightarrow 5.4 \pm 0.1
164			\mathrm{ms}^{-1}$ for the secondary prism, we found no significant change in
165			the interfacial width. This follows our previous findings of the basal and
166			prismatic systems, in which the interfacial width was invarient of the
167			shear rate of the ice. The interfacial width of the quiescent basal and
168			prismatic systems was found to be $3.2 \pm 0.4$ \AA\ and $3.6 \pm 0.2$ \AA\
169			respectively. Over the range of shear rates investigated, $0.6 \pm 0.3
170			\mathrm{ms}^{-1} \rightarrow 5.3 \pm 0.5 \mathrm{ms}^{-1}$ for the basal
171			system and $0.9 \pm 0.2 \mathrm{ms}^{-1} \rightarrow 4.5 \pm 0.1
172	plouden	4194	\mathrm{ms}^{-1}$ for the prismatic.
173
174			These results indicate that the surface structure of the exposed ice crystal
175			has little to no effect on how far into the bulk the ice-like structural
176			ordering is. Also, it appears that the interface is not structurally effected
177			by shearing the ice through water.
178
179
180	plouden	4192	\subsection{Orientational dynamics}
181	plouden	4215	%Should we include the math here?
182			The orientational time correlation function,
183			\begin{equation}\label{C(t)1}
184			C_{2}(t)=\langle P_{2}(\mathbf{u}(0)\cdot \mathbf{u}(t))\rangle,
185			\end{equation}
186			helps indicate the local environment around the water molecules. The function
187			begins with an initial value of unity, and decays to zero as the water molecule
188			loses memory of its former orientation. Observing the rate at which this decay
189			occurs can provide insight to the mechanism and timescales for the relaxation.
190			In eq. \eqref{C(t)1}, $P_{2}$ is the second-order Legendre polynomial, and
191			$\mathbf{u}$ is the the bisecting HOH vector. The angle brackets indicate
192			an ensemble average over all the water molecules in a given spatial region.
193
194	plouden	4194	To investigate the dynamics of the water molecules across the interface, the
195	plouden	4215	systems were divided in the $z$-dimension into bins, each $\approx$ 3 \AA\
196			wide, and \eqref{C(t)1} was computed for each of the bins. A water
197			molecule was allocated to a particular bin if it was initially in the bin
198			at time zero. To compute \eqref{C(t)1}, each 0.5 ns simulation was followed
199			by an additional 200 ps microcanonical (NVE) simulation during which the
200			position and orientations of each molecule were recorded every 0.1 ps.
201
202			The data obtained for each bin was then fit to a triexponential decay given by
203			\begin{equation}\label{C(t)_fit}
204			C_{2}(t) \approx a e^{-t/\tau_\mathrm{short}} + b e^{-t/\tau_\mathrm{middle}} +\
205			c
206			e^{-t/\tau_\mathrm{long}} + (1-a-b-c)
207			\end{equation}
208			where $\tau_{short}$ corresponds to the librational motion of the water
209			molecules, $\tau_{middle}$ corresponds to jumps between the breaking and
210			making of hydrogen bonds, and $\tau_{long}$ corresponds to the translational
211			motion of the water molecules. The last term in \eqref{C(t)_fit} accounts
212			for the water molecules trapped in the ice which do not experience any
213			long-time orientational decay.
214	plouden	4192
215	plouden	4215	In Figures \ref{fig:PyrOrient} and \ref{fig:SPorient} we see the $z$-coordinate
216			profiles for the three decay constants, $\tau_{short}$ (panel a),
217			$\tau_{middle}$ (panel b),
218			and $\tau_{long}$ (panel c) for the pyramidal and secondary prismatic systems
219			respectively. The control experiments (no shear) are shown in black, and
220			an experiment with an imposed momentum flux is shown in red. The vertical
221			dotted line traversing all three panels denotes the midpoint of the
222			interface as determined by the local tetrahedral order parameter fitting.
223			In the liquid regions of both systems, we see that $\tau_{middle}$ and
224			$\tau_{long}$ have approximately consistent values of $3-6$ ps and $30-40$ ps,
225			resepctively, and increase in value as we approach the interface. Conversely,
226			in panel a, we see that $\tau_{short}$ decreases from the liquid value
227			of $72-76$ fs as we approach the interface. We believe this speed up is due to
228			the constrained motion of librations closer to the interface. Both the
229			approximate values for the decays and relative trends match those reported
230			previously for the basal and prismatic interfaces.
231	plouden	4192
232	plouden	4215	As done previously, we have attempted to quantify the distance, $d_{pyramidal}$
233			and $d_{secondary prism}$, from the
234			interface that the deviations from the bulk liquid values begin. This was done
235			by fitting the orientational decay constant $z$-profiles by
236			\begin{equation}\label{tauFit}
237			\tau(z)\approx\tau_{liquid}+(\tau_{solid}-\tau_{liquid})e^{-(z-z_{wall})/d}
238			\end{equation}
239			where $\tau_{liquid}$ and $\tau_{solid}$ are the liquid and projected solid
240			values of the decay constants, $z_{wall}$ is the location of the interface,
241			and $d$ is the displacement from the interface at which these deviations
242			occur. The values for $d_{pyramidal}$ and $d_{secondary prismatic}$ were
243			determined
244			for each of the decay constants, and then averaged for better statistics
245			($\tau_{middle}$ was ommitted for secondary prism). For the pyramidal system,
246			$d_{pyramidal}$ was found to be 2.7 \AA\ for both the control and the sheared
247			system. We found $d_{secondary prismatic}$ to be slightly larger than
248			$d_{pyramidal}$ for both the control and with an applied shear, with
249			displacements of $4$ \AA\ for the control system and $3$ \AA\ for the
250			experiment with the imposed momentum flux. These values are consistent with
251			those found for the basal ($d_{basal}\approx2.9$ \AA\ ) and prismatic
252			($d_{prismatic}\approx3.5$ \AA\ ) systems.
253	plouden	4192
254	plouden	4194	\subsection{Coefficient of friction of the interfaces}
255	plouden	4215	While investigating the kinetic coefficient of friction for the larger
256			prismatic system, there was found to be a dependence for $\mu_k$
257			on the temperature of the liquid water in the system. We believe this
258			dependence
259			arrises from the sharp discontinuity of the viscosity for the SPC/E model
260			at temperatures approaching 200 K\cite{kuang12}. Due to this, we propose
261			a weighting to the structural interfacial parameter, $\kappa$ by the
262			viscosity at $225$ K, the temperature of the interface. $\kappa$ is
263			traditionally defined as
264			\begin{equation}\label{kappa}
265			\kappa = \eta/\delta
266			\end{equation}
267			where $\eta$ is the viscosity and $\delta$ is the slip length.
268			In our ice/water shearing simulations, the system has reached a steady state
269			when the applied force,
270	plouden	4194
271	plouden	4215	\begin{equation}
272			f_{applied} = \mathbf{j}_z(\mathbf{p})L_x L_y
273			\end{equation}
274			is equal to the frictional force resisting the motion of the ice block
275			\begin{equation}
276			f_{friction} = \frac{\eta \mathbf{v} L_x L_y}{\delta}
277			\end{equation}
278			where $\mathbf{v}$ is the relative velocity of the liquid from the ice.
279			When this condition is met, we are able to solve the resulting expression to
280			obtain,
281			\begin{equation}\label{force_equality}
282			\frac{\eta}{\delta} = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}}
283			\end{equation}
284			From \eqref{kappa}, \eqref{force_equality} becomes
285			\begin{equation}
286			\kappa = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}}
287			\end{equation}
288			which we will multiply by a viscosity weighting term to reach
289			\begin{equation} \label{kappa2}
290			\kappa = \frac{\mathbf{j}_z(\mathbf{p})}{\mathbf{v}} \frac{\eta(225)}{\eta(T)}
291			\end{equation}
292			Assuming linear response theory is valid, an expression for ($\eta$) can
293			be found from the imposed momentum flux and the measured velocity gradient.
294			\begin{equation}\label{eta_eq}
295			\eta = \frac{\mathbf{j}_z(\mathbf{p})}{\frac{\partial v_x}{\partial z}}
296			\end{equation}
297			Substituting eq \eqref{eta_eq} into eq \eqref{kappa2} we arrive at
298			\begin{equation}
299			\kappa = \frac{\frac{\partial v_x}{\partial z}}{\mathbf{v}}\eta(225)
300			\end{equation}
301	plouden	4194
302	plouden	4219	\begin{table}[h]
303			\centering
304			\caption{$\kappa$ values for the basal, prismatic, pyramidal, and secondary prismatic facets of Ice-I$_\mathrm{h}$}
305			\label{tab:kappa}
306			\begin{tabular}{\|ccc\|} \hline
307			& \multicolumn{2}{c\|}{$\kappa_{Drag direction}$} \\
308			Interface & $\kappa_{x}$ & $\kappa_{y}$ \\ \hline
309			basal & $0.00059 \pm 0.00003$ & $0.00065 \pm 0.00008$ \\
310			prismatic & $0.00030 \pm 0.00002$ & $0.00030 \pm 0.00001$ \\
311			pyramidal & $0.00058 \pm 0.00004$ & $0.00061 \pm 0.00005$ \\
312			secondary prism & $0.00035 \pm 0.00001$ & $0.00033 \pm 0.00002$ \\ \hline
313			\end{tabular}
314			\end{table}
315
316
317	plouden	4215	To obtain the value of $\eta(225)$ for the SPC/E model, a $31.09 \times 29.38
318			\times 124.39$ \AA\ box with 3744 water molecules was equilibrated to 225K,
319	plouden	4219	and 5 unique shearing experiments were performed. Each experiment was
320			conducted in the microcanonical ensemble (NVE) and were 5 ns in
321			length. The VSS were attempted every timestep, which was set to 2 fs.
322			For our SPC/E systems, we found $\eta(225)$ to be $0.0148 \pm 0.0007 Pa s$,
323			roughly ten times larger than the value found for 280 K SPC/E water by
324			Kuang\cite{kuang12}.
325	plouden	4215
326
327	plouden	4194	\begin{table}[h]
328			\centering
329			\caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript\
330	gezelter	4217	{-1}). \\
331			\textsuperscript{a} See ref. \onlinecite{Louden13}. }
332	plouden	4194	\label{tab:lambda}
333			\begin{tabular}{\|ccc\|} \hline
334			& \multicolumn{2}{c\|}{Drag direction} \\
335			Interface & $x$ & $y$ \\ \hline
336			basal\textsuperscript{a} & $0.08 \pm 0.02$ & $0.09 \pm 0.03$ \\
337	plouden	4215	prismatic (T = 225)\textsuperscript{a} & $0.037 \pm 0.008$ & $0.04 \pm 0.01$ \\
338			prismatic (T = 230) & $0.10 \pm 0.01$ & $0.070 \pm 0.006$\\
339	plouden	4194	pyramidal & $0.13 \pm 0.03$ & $0.14 \pm 0.03$ \\
340			secondary prism & $0.13 \pm 0.02$ & $0.12 \pm 0.03$ \\ \hline
341			\end{tabular}
342			\end{table}
343
344
345	plouden	4192	\begin{figure}
346			\includegraphics[width=\linewidth]{Pyr-orient}
347			\caption{\label{fig:PyrOrient} The three decay constants of the
348			orientational time correlation function, $C_2(t)$, for water as a function
349			of distance from the center of the ice slab. The vertical dashed line
350			indicates the edge of the pyramidal ice slab determined by the local order
351			tetrahedral parameter. The control (black circles) and sheared (red squares)
352			experiments were fit by a shifted exponential decay (Eq. 9\cite{Louden13})
353			shown by the black and red lines respectively. The upper two panels show that
354			translational and hydrogen bond making and breaking events slow down
355			through the interface while approaching the ice slab. The bottom most panel
356			shows the librational motion of the water molecules speeding up approaching
357			the ice block due to the confined region of space allowed for the molecules
358			to move in.}
359			\end{figure}
360
361			\begin{figure}
362			\includegraphics[width=\linewidth]{SP-orient-less}
363			\caption{\label{fig:SPorient} Decay constants for $C_2(t)$ at the secondary
364			prism face. Panel descriptions match those in \ref{fig:PyrOrient}.}
365			\end{figure}
366
367
368
369			\section{Conclusion}
370			Conclude conclude conclude...
371
372
373	gezelter	4217	\begin{acknowledgments}
374			Support for this project was provided by the National
375			Science Foundation under grant CHE-1362211. Computational time was
376			provided by the Center for Research Computing (CRC) at the
377			University of Notre Dame.
378			\end{acknowledgments}
379	plouden	4192
380			\newpage
381	gezelter	4217
382	plouden	4192	\bibliography{iceWater}
383
384			\end{document}