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}
In this follow up paper of the basal and prismatic facets of the 
Ice-I$_\mathrm{h}$/water interface, we present the 
pyramidal and secondary prismatic 
interfaces for both the quiescent and sheared systems. The structural and
dynamic interfacial widths for all four crystal facets were found to be in good
agreement, and were found to be independent of the shear rate over the shear 
rates investigated. 
Decomposition of the molecular orientational time correlation function showed
different behavior for the short- and longer-time decay components approaching
normal to the interface. Lastly we show through calculation of the interfacial
friction coefficient that the basal and pyramidal facets are more 
hydrophilic than the prismatic and secondary prismatic facets.

\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 our 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 used by Bryk and Haymet in a previous study
of ice/water interfaces.\cite{Bryk2004b} 

The local tetrahedral order parameter, $q(z)$, is given by
\begin{equation}
q(z) = \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3}
\sum_{j=i+1}^{4} \bigg(\cos\psi_{ikj}+\frac{1}{3}\bigg)^2\Bigg)
\delta(z_{k}-z)\mathrm{d}z \Bigg/ N_z
\label{eq:qz}
\end{equation}
where $\psi_{ikj}$ is the angle formed between the oxygen sites of molecules
$i$,$k$, and $j$, where the centeral oxygen is located within molecule $k$ and 
molecules $i$ and $j$ are two of the closest four water molecules
around molecule $k$. All four closest neighbors of molecule $k$ are also
required to reside within the first peak of the pair distribution function
for molecule $k$ (typically $<$ 3.41 \AA\ for water). 
$N_z = \int\delta(z_k - z) \mathrm{d}z$ is a normalization factor to account 
for the varying population of molecules within each finite-width bin.

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 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, 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 interfacial friction coefficient, $\kappa$ by the 
shear viscosity at 225 K. The interfacial friction coefficient relates
the shear stress with the relative velocity of the fluid normal to the 
interface:
\begin{equation}\label{Shenyu-13}
j_{z}(p_{x}) = \kappa[v_{x}(fluid)-v_{x}(solid)]
\end{equation}
where $j_{z}(p_{x})$ is the applied momentum flux (shear stress) across $z$ 
in the 
$x$-dimension, and $v_{x}$(fluid) and $v_{x}$(solid) are the velocities
directly adjacent to the interface. The shear viscosity, $\eta(T)$, of the 
fluid can be determined if we assume a linear response of the momentum 
gradient to the applied shear stress by
\begin{equation}\label{Shenyu-11}
j_{z}(p_{x}) = \eta(T) \frac{\partial v_{x}}{\partial z}.
\end{equation}
Using eqs \eqref{Shenyu-13} and \eqref{Shenyu-11}, we can find the following
expression for $\kappa$,
\begin{equation}\label{kappa-1}
\kappa = \eta(T) \frac{\partial v_{x}}{\partial z}\frac{1}{[v_{x}(fluid)-v_{x}(solid)]}.
\end{equation}
Here is where we will introduce the weighting term of $\eta(225)/\eta(T)$ 
giving us
\begin{equation}\label{kappa-2}
\kappa = \frac{\eta(225)}{[v_{x}(fluid)-v_{x}(solid)]}\frac{\partial v_{x}}{\partial z}.
\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 SPC/E 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 bulk water by
Kuang\cite{kuang12}. 

The interfacial friction coefficient, $\kappa$, can equivalently be expressed 
as the ratio of the viscosity of the fluid to the slip length, $\delta$, which
is an indication of how 'slippery' the interface is.
\begin{equation}\label{kappa-3}
\kappa = \frac{\eta}{\delta}
\end{equation} 
In each of the systems, the interfacial temperature was kept fixed to 225K, 
which ensured the viscosity of the fluid at the
interace was approximately the same. Thus, any significant variation in
$\kappa$ between
the systems indicates differences in the 'slipperiness' of the interfaces.
As each of the ice systems are sheared relative to liquid water, the
'slipperiness' of the interface can be taken as an indication of how
hydrophobic or hydrophilic the interface is. The calculated $\kappa$ values
found for the four crystal facets of Ice-I$_\mathrm{h}$ investigated are shown
in Table \ref{tab:kapa}. The basal and pyramidal facets were found to have 
similar values of $\kappa \approx$ 0.0006 
(amu \AA\ \textsuperscript{-2} fs\textsuperscript{-1}), while $\kappa \approx$ 
0.0003 (amu \AA\ \textsuperscript{-2} fs\textsuperscript{-1}) were found for 
the prismatic and secondary prismatic systems.
These results indicate that the prismatic and secondary prismatic facets are
more hydrophobic than the basal and pyramidal facets.
%This indicates something about the similarity between the two facets that 
%share similar values...
%Maybe find values for kappa for other materials to compare against?

\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}$ (amu \AA\ \textsuperscript{-2} fs\textsuperscript{-1})} \\
 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}


%\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}
We present the results of molecular dynamics simulations of the pyrmaidal
and secondary prismatic facets of an SPC/E model of the 
Ice-I$_\mathrm{h}$/water interface. The ice was sheared through the liquid
water while being exposed to a thermal gradient to maintain a stable 
interface by using the minimal perturbing VSS RNEMD method. In agreement with 
our previous findings for the basal and prismatic facets, the interfacial 
width was found to be independent of shear rate as measured by the local 
tetrahedral order parameter. This width was found to be 
3.2~$\pm$0.2~\AA\ for both the pyramidal and the secondary prismatic systems.
These values are in good agreement with our previously calculated interfacial
widths for the basal (3.2~$\pm$0.4~\AA\ ) and prismatic (3.6~$\pm$0.2~\AA\ )
systems. 

Orientational dynamics of the Ice-I$_\mathrm{h}$/water interfaces were studied
by calculation of the orientational time correlation function at varying
displacements normal to the interface. The decays were fit
to a tri-exponential decay, where the three decay constants correspond to
the librational motion of the molecules driven by the restoring forces of 
existing hydrogen bonds ($\tau_{short}$ $\mathcal{O}$(10 fs)), jumps between 
two different hydrogen bonds ($\tau_{middle}$ $\mathcal{O}$(1 ps)), and 
translational motion of the molecules ($\tau_{long}$ $\mathcal{O}$(100 ps)).
$\tau_{short}$ was found to decrease approaching the interface due to the 
constrained motion of the molecules as the local environment becomes more
ice-like. Conversely, the two longer-time decay constants were found to 
increase at small displacements from the interface. As seen in our previous 
work on the basal and prismatic facets, there appears to be a dynamic
interface width at which deviations from the bulk liquid values occur.
We had previously found $d_{basal}$ and $d_{prismatic}$ to be approximately 
2.8~\AA\ and 3.5~\AA\~. We found good agreement of these values for the 
pyramidal and secondary prismatic systems with $d_{pyramidal}$ and
$d_{secondary prism}$ to be 2.7~\AA\ and 3~\AA\~. For all of the facets, there
was found to be no apparent dependence of the dynamic width on the shear rate.
  
%Paragraph summarizing the Kappa values
The interfacial friction coefficient, $\kappa$, was determined for each of the
interfaces. We were able to reach an expression for $\kappa$ as a function of
the velocity profile of the system and is scaled by the viscosity of the liquid
at 225 K. In doing so, we have obtained an expression for $\kappa$ which is 
independent of temperature differences of the liquid water at far displacements
from the interface. We found the basal and pyramidal facets to have
similar $\kappa$ values, of $\kappa \approx$ 0.0006 
(amu \AA\ \textsuperscript{-2} fs\textsuperscript{-1}). However, the
prismatic and secondary prismatic facets were found to have $\kappa$ values of 
$\kappa \approx$ 0.0003 (amu \AA\ \textsuperscript{-2} fs\textsuperscript{-1}).
As these ice facets are being sheared relative to liquid water, with the
structural and dynamic width of all four 
interfaces being approximately the same, the difference in the coefficient of 
friction indicates the hydrophilicity of the crystal facets are not 
equivalent. Namely, that the basal and pyramidal facets of Ice-I$_\mathrm{h}$
are more hydrophilic than the prismatic and secondary prismatic facets.


\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:	4225
Committed:	Tue Sep 16 20:27:18 2014 UTC (9 years, 11 months ago) by plouden
Content type:	application/x-tex
File size:	24102 byte(s)
Log Message:	Wrote an Abstract