trunk/iceWater2/iceWater2.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
%\usepackage{endfloat}                                                        
\usepackage{caption}
%\usepackage{epsf}                                                            
%\usepackage{tabularx}                                                        
\usepackage{graphicx}
\usepackage{multirow}
\usepackage{wrapfig}
%\usepackage{booktabs}                                                        
%\usepackage{bibentry}                                                        
%\usepackage{mathrsfs}                                                        
%\usepackage[ref]{overcite}                                                   
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures                                     
%\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}                 
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

%\renewcommand\citemid{\ } % no comma in optional referenc note               
\bibpunct{}{}{,}{s}{}{;}
\bibliographystyle{aip}


% \documentclass[journal = jpccck, manuscript = article]{achemso}            
% \setkeys{acs}{usetitle = true}                                              
% \usepackage{achemso}                                                        
% \usepackage{natbib}                                                         
% \usepackage{multirow}                                                       
% \usepackage{wrapfig}                                                        
% \usepackage{fixltx2e}                                                        

\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions                                                                  
\usepackage{url}


\begin{document}

\title{Simulations of solid-liquid friction at Secondary Prism and Pyramidal ice-I$_\mathrm{h}$ / water interfaces}

\author{Patrick B. Louden and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail:
  gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}
\maketitle
\begin{doublespace}

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

\newpage

\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 structured differently 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,Hayword01,Hayword02,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} 

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 6 unique shearing experiments were performed. Each experiment was
conducted in the microcanonical ensemble (NVE) and were 1 ns in 
length. The VSS were attempted every timestep, which was set to 2 fs.  


\begin{table}[h]
\centering
\caption{Solid-liquid friction coefficients (measured in amu~fs\textsuperscript\
{-1}) }
\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}
\caption{\textsuperscript{a}Reference \cite{Louden13}}
\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...

\section{Acknowledgements}
Support for this progect was provided by the National Science Foundation under grant CHE-0848243. Computational time was provided by the Center for Research Computing (CRC) at the University of Notre Dame.


\newpage
\bibliography{iceWater}

\end{doublespace}

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