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