| 11 |
|
\usepackage{url} |
| 12 |
|
|
| 13 |
|
|
| 14 |
< |
\title{Do the facets of ice $I_\mathrm{h}$ crystals have different |
| 15 |
< |
friction coefficients? Simulating shear in ice/water interfaces} |
| 14 |
> |
\title{Solid-liquid friction at ice-I$_\mathrm{h}$ / water interfaces} |
| 15 |
|
|
| 16 |
|
\author{P. B. Louden} |
| 17 |
|
\author{J. Daniel Gezelter} |
| 25 |
|
\begin{document} |
| 26 |
|
|
| 27 |
|
\begin{abstract} |
| 28 |
< |
We have investigated the structural properties of the basal and |
| 29 |
< |
prismatic facets of an SPC/E model of the ice Ih / water interface |
| 30 |
< |
when the solid phase is being drawn through liquid water (i.e. sheared |
| 31 |
< |
relative to the fluid phase). To impose the shear, we utilized a |
| 32 |
< |
reverse non-equilibrium molecular dynamics (RNEMD) method that creates |
| 33 |
< |
non-equilibrium conditions using velocity shearing and scaling (VSS) |
| 34 |
< |
moves of the molecules in two physically separated slabs in the |
| 35 |
< |
simulation cell. This method can create simultaneous temperature and |
| 36 |
< |
velocity gradients and allow the measurement of friction transport |
| 37 |
< |
properties at interfaces. We found the interfacial width for the basal and prismatic facets to be 3.2$\pm$0.4 \AA\ and 3.6$\pm$0.2 \AA\ with no momentum gradient present. The interfacial width was found to be independent of shear rate over the range of shear rates investigated, 0.6$\pm$0.3 ms\textsuperscript{-1} to 5.3$\pm$0.5 ms\textsuperscript{-1} for the basal system and 0.9$\pm$0.2 ms\textsuperscript{-1} to 4.5$\pm$0.1 ms\textsuperscript{-1} for the prismatic. The orientational time correlation function was evaluated at varying displacements from the interface, and fit to a triexponential decay with decay times $\tau_{short}$, $\tau_{middle}$, and $\tau_{long}$. We found that each decay time was independent of the shear rate, and that they each deviated from a bulk liquid value at the same displacement $d_{basal}$ and $d_{prismatic}$. We determined $d_{basal}$ and $d_{prismatic}$ to be approximately 15 \AA\ and 8 \AA\ . Both $\tau_{middle}$ and $\tau_{long}$ increased in value in decreasing displacements from $d_{basal}$ and $d_{prismatic}$, while we found that $\tau_{short}$ had the inverse correlation. Lastly we present calculations of the friction coefficients ($\lambda$) for the basal and prismatic facets. We have found $\lambda_{basal}$ to be (0.07$\pm$0.01) amu fs\textsuperscript{-1} and $\lambda_{prismatic}$ to be (0.032$\pm$0.007) amu fs\textsuperscript{-1}. |
| 28 |
> |
We have investigated the structural and dynamic properties of the |
| 29 |
> |
basal and prismatic facets of an ice I$_\mathrm{h}$ / water |
| 30 |
> |
interface when the solid phase is being drawn through the liquid |
| 31 |
> |
(i.e. sheared relative to the fluid phase). To impose the shear, we |
| 32 |
> |
utilized a velocity-shearing and scaling (VSS) approach to reverse |
| 33 |
> |
non-equilibrium molecular dynamics (RNEMD). This method can create |
| 34 |
> |
simultaneous temperature and velocity gradients and allow the |
| 35 |
> |
measurement of transport properties at interfaces. The interfacial |
| 36 |
> |
width was found to be independent of relative velocity of the ice |
| 37 |
> |
and liquid layers over a wide range of shear rates. Decays of |
| 38 |
> |
molecular orientational time correlation functions for gave very |
| 39 |
> |
similar estimates for the width of the interfaces, although the |
| 40 |
> |
short- and longer-time decay components of the orientational |
| 41 |
> |
correlation functions behave differently closer to the interface. |
| 42 |
> |
Although both facets of ice are in ``stick'' boundary conditions in |
| 43 |
> |
liquid water, the solid-liquid friction coefficient was found to be |
| 44 |
> |
different for the basal and prismatic facets of ice. |
| 45 |
|
\end{abstract} |
| 46 |
|
|
| 47 |
|
\newpage |
| 56 |
|
|
| 57 |
|
%Gay02: cites many other ice/water papers, make sure to cite them. |
| 58 |
|
|
| 59 |
< |
Understanding the ice/water interface is essential for explaining complex processes such as nucleartion and crystal growth\cite{Han92,Granasy95,Vanfleet95}, crystal melting\cite{Weber83,Han92,Sakai96,Sakai96B}, and biological interfacial processes, such as the antifreeze protein found in winter flounder\cite{Wierzbicki07, Chapsky97}. These processes have been studied at the fundamental level of the ice/water interface by several groups, including studying the structure and width of the interface. Haymet \emph{et al.} have done extensive work on ice Ih, the most common form of ice on Earth, including characterizing and determining the width of the ice/water interface for the SPC, SPC/E, CF1, and TIP4P models for water. \cite{Karim87,Karim88,Karim90,Hayward01,Bryk02,Hayward02,Gay02} More recently, Haymet \emph{et al.} have been investigating the effects cations and anions have on crystal nucleation\cite{Bryk04,Smith05,Wilson08,Wilson10}. Nada \emph{et al.} have also studied the ice/water interface\cite{Nada95,Nada00,Nada03,Nada12}. They have found that the different facets of ice Ih have different growth rates, primarily, that the prismatic facet grows faster than the basal facet due to the mechanism of the crystal growth being the reordering of the hydrogen bonding network\cite{Nada05}. |
| 59 |
> |
Understanding the ice/water interface is essential for explaining |
| 60 |
> |
complex processes such as nucleation and crystal |
| 61 |
> |
growth,\cite{Han92,Granasy95,Vanfleet95} crystal |
| 62 |
> |
melting,\cite{Weber83,Han92,Sakai96,Sakai96B} and some fascinating |
| 63 |
> |
biological processes, such as the behavior of the antifreeze proteins |
| 64 |
> |
found in winter flounder,\cite{Wierzbicki07, Chapsky97} and certain |
| 65 |
> |
terrestrial arthropods.\cite{Duman:2001qy,Meister29012013} There has |
| 66 |
> |
been significant progress on understanding the structure and dynamics |
| 67 |
> |
of quiescent ice/water interfaces utilizing both theory and |
| 68 |
> |
experiment. Haymet \emph{et al.} have done extensive work on ice I$_\mathrm{h}$, |
| 69 |
> |
including characterizing and determining the width of the ice/water |
| 70 |
> |
interface for the SPC, SPC/E, CF1, and TIP4P models for |
| 71 |
> |
water.\cite{Karim87,Karim88,Karim90,Hayward01,Bryk02,Hayward02,Gay02} |
| 72 |
> |
More recently, Haymet \emph{et al.} have investigated the effects |
| 73 |
> |
cations and anions have on crystal |
| 74 |
> |
nucleation\cite{Bryk04,Smith05,Wilson08,Wilson10}. Nada \emph{et al.} |
| 75 |
> |
have also studied ice/water |
| 76 |
> |
interfaces,\cite{Nada95,Nada00,Nada03,Nada12} and have found that the |
| 77 |
> |
differential growth rates of the facets of ice I$_\mathrm{h}$ are due to the the |
| 78 |
> |
reordering of the hydrogen bonding network\cite{Nada05}. |
| 79 |
|
|
| 80 |
< |
Another complex process which requires investigation at the ice/water interface is the movement of water over ice, such as icebergs floating in the ocean. In addition to understanding the structure and width of the interface, it is pertinent to understand the friction caused by the shearing of water across the ice to understand this process. However, until recently, simulations of this nature were not possible. |
| 81 |
< |
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. In this paper, we investigate the width and the friction coefficient of the ice/water interface as the ice is sheared through the liquid under a weak thermal gradient. |
| 80 |
> |
The movement of liquid water over the facets of ice has been less |
| 81 |
> |
thoroughly studied than the quiescent surfaces. This process is |
| 82 |
> |
potentially important in understanding transport of large blocks of |
| 83 |
> |
ice in water (which has important implications in the earth sciences), |
| 84 |
> |
as well as the relative motion of crystal-crystal interfaces that have |
| 85 |
> |
been separated by nanometer-scale fluid domains. In addition to |
| 86 |
> |
understanding both the structure and thickness of the interfacial |
| 87 |
> |
regions, it is important to understand the molecular origin of |
| 88 |
> |
friction, drag, and other changes in dynamical properties of the |
| 89 |
> |
liquid in the regions close to the surface that are altered by the |
| 90 |
> |
presence of a shearing of the bulk fluid relative to the solid phase. |
| 91 |
|
|
| 92 |
+ |
In this work, we apply a recently-developed velocity shearing and |
| 93 |
+ |
scaling approach to reverse non-equilibrium molecular dynamics |
| 94 |
+ |
(VSS-RNEMD). This method makes it possible to calculate transport |
| 95 |
+ |
properties like the interfacial thermal conductance across |
| 96 |
+ |
heterogeneous interfaces,\cite{Kuang12} and can create simultaneous |
| 97 |
+ |
temperature and velocity gradients and allow the measurement of |
| 98 |
+ |
friction and thermal transport properties at interfaces. This has |
| 99 |
+ |
allowed us to investigate the width of the ice/water interface as the |
| 100 |
+ |
ice is sheared through the liquid, while simultaneously imposing a |
| 101 |
+ |
weak thermal gradient to prevent frictional heating of the interface. |
| 102 |
+ |
In the sections that follow, we discuss the methodology for creating |
| 103 |
+ |
and simulating ice/water interfaces under shear and provide results |
| 104 |
+ |
from both structural and dynamical correlation functions. We also |
| 105 |
+ |
show that the solid-liquid interfacial friction coefficient depends |
| 106 |
+ |
sensitively on the details of the surface morphology. |
| 107 |
+ |
|
| 108 |
|
\section{Methodology} |
| 109 |
|
|
| 110 |
< |
\subsection{Stable ice I$_\mathrm{h}$ / water interfaces} |
| 110 |
> |
\subsection{Stable ice I$_\mathrm{h}$ / water interfaces under shear} |
| 111 |
|
|
| 112 |
< |
The structure of ice I$_\mathrm{h}$ is well understood; it |
| 112 |
> |
The structure of ice I$_\mathrm{h}$ is very well understood; it |
| 113 |
|
crystallizes in a hexagonal space group P$6_3/mmc$, and the hexagonal |
| 114 |
|
crystals of ice have two faces that are commonly exposed, the basal |
| 115 |
|
face $\{0~0~0~1\}$, which forms the top and bottom of each hexagonal |
| 116 |
|
plate, and the prismatic $\{1~0~\bar{1}~0\}$ face which forms the |
| 117 |
|
sides of the plate. Other less-common, but still important, faces of |
| 118 |
< |
ice I$_\mathrm{h}$ are the secondary prism face, $\{1~1~\bar{2}~0\}$, |
| 119 |
< |
and the prismatic face, $\{2~0~\bar{2}~1\}$. |
| 118 |
> |
ice I$_\mathrm{h}$ are the secondary prism, $\{1~1~\bar{2}~0\}$, and |
| 119 |
> |
the prismatic, $\{2~0~\bar{2}~1\}$, faces. Ice I$_\mathrm{h}$ is |
| 120 |
> |
normally proton disordered in bulk crystals, although the surfaces |
| 121 |
> |
probably have a preference for proton ordering along strips of |
| 122 |
> |
dangling H-atoms and Oxygen lone pairs.\cite{Buch:2008fk} |
| 123 |
|
|
| 71 |
– |
Ice I$_\mathrm{h}$ is normally proton disordered in bulk crystals, |
| 72 |
– |
although the surfaces probably have a preference for proton ordering |
| 73 |
– |
along strips of dangling H-atoms and Oxygen lone |
| 74 |
– |
pairs.\cite{Buch:2008fk} |
| 75 |
– |
|
| 124 |
|
For small simulated ice interfaces, it is useful to have a |
| 125 |
|
proton-ordered, but zero-dipole crystal that exposes these strips of |
| 126 |
< |
dangling H-atoms and lone pairs. Also, if we're going to place |
| 127 |
< |
another material in contact with one of the ice crystalline planes, it |
| 128 |
< |
is useful to have an orthorhombic (rectangular) box to work with. A |
| 129 |
< |
recent paper by Hirsch and Ojam\"{a}e describes how to create |
| 130 |
< |
proton-ordered bulk ice I$_\mathrm{h}$ in alternative orthorhombic |
| 83 |
< |
cells.\cite{Hirsch04} |
| 126 |
> |
dangling H-atoms and lone pairs. When placing another material in |
| 127 |
> |
contact with one of the ice crystalline planes, it is useful to have |
| 128 |
> |
an orthorhombic (rectangular) box. A recent paper by Hirsch and |
| 129 |
> |
Ojam\"{a}e describes how to create proton-ordered bulk ice |
| 130 |
> |
I$_\mathrm{h}$ in alternative orthorhombic cells.\cite{Hirsch04} |
| 131 |
|
|
| 132 |
|
We are using Hirsch and Ojam\"{a}e's structure 6 which is an |
| 133 |
|
orthorhombic cell ($P2_12_12_1$) that produces a proton-ordered |
| 134 |
< |
version of ice Ih. Table \ref{tab:equiv} contains a mapping between |
| 134 |
> |
version of ice I$_\mathrm{h}$. Table \ref{tab:equiv} contains a mapping between |
| 135 |
|
the Miller indices in the P$6_3/mmc$ crystal system and those in the |
| 136 |
|
Hirsch and Ojam\"{a}e $P2_12_12_1$ system. |
| 137 |
|
|
| 138 |
|
\begin{wraptable}{r}{3.5in} |
| 139 |
+ |
\caption{Mapping between the Miller indices of four facets of ice in |
| 140 |
+ |
the $P6_3/mmc$ crystal system to the orthorhombic $P2_12_12_1$ |
| 141 |
+ |
system in reference \protect\cite{Hirsch04}} |
| 142 |
+ |
\label{tab:equiv} |
| 143 |
|
\begin{tabular}{|ccc|} \hline |
| 144 |
|
& hexagonal & orthorhombic \\ |
| 145 |
|
& ($P6_3/mmc$) & ($P2_12_12_1$) \\ |
| 153 |
|
|
| 154 |
|
Structure 6 from the Hirsch and Ojam\"{a}e paper has lattice |
| 155 |
|
parameters $a = 4.49225$, $b = 7.78080$, $c = 7.33581$ and two water |
| 156 |
< |
molecules whose atoms reside at the following fractional coordinates: |
| 156 |
> |
molecules whose atoms reside at fractional coordinates given in table |
| 157 |
> |
\ref{tab:p212121}. To construct the basal and prismatic interfaces, |
| 158 |
> |
these crystallographic coordinates were used to construct an |
| 159 |
> |
orthorhombic unit cell which was then replicated in all three |
| 160 |
> |
dimensions yielding a proton-ordered block of ice I$_{h}$. To expose |
| 161 |
> |
the desired face, the orthorhombic representation was then cut along |
| 162 |
> |
the ($001$) or ($100$) planes for the basal and prismatic faces |
| 163 |
> |
respectively. The resulting block was rotated so that the exposed |
| 164 |
> |
faces were aligned with the $z$-dimension normal to the exposed face. |
| 165 |
> |
The block was then cut along two perpendicular directions in a way |
| 166 |
> |
that allowed for perfect periodic replication in the $x$ and $y$ axes, |
| 167 |
> |
creating a slab with either the basal or prismatic faces exposed along |
| 168 |
> |
the $z$ axis. The slab was then replicated in the $x$ and $y$ |
| 169 |
> |
dimensions until a desired sample size was obtained. |
| 170 |
|
|
| 171 |
|
\begin{wraptable}{r}{3.25in} |
| 172 |
< |
\begin{tabular}{|ccccc|} \hline |
| 173 |
< |
atom label & type & x & y & z \\ \hline |
| 174 |
< |
O$_{a}$ & O & 0.75 & 0.1667 & 0.4375 \\ |
| 175 |
< |
H$_{1a}$ & H & 0.5735 & 0.2202 & 0.4836 \\ |
| 176 |
< |
H$_{2a}$ & H & 0.7420 & 0.0517 & 0.4836 \\ |
| 177 |
< |
O$_{b}$ & O & 0.25 & 0.6667 & 0.4375 \\ |
| 178 |
< |
H$_{1b}$ & H & 0.2580 & 0.6693 & 0.3071 \\ |
| 179 |
< |
H$_{2b}$ & H & 0.4265 & 0.7255 & 0.4756 \\ \hline |
| 172 |
> |
\caption{Fractional coordinates for water in the orthorhombic |
| 173 |
> |
$P2_12_12_1$ system for ice I$_\mathrm{h}$ in reference |
| 174 |
> |
\protect\cite{Hirsch04}} |
| 175 |
> |
\label{tab:p212121} |
| 176 |
> |
\begin{tabular}{|cccc|} \hline |
| 177 |
> |
atom type & x & y & z \\ \hline |
| 178 |
> |
O & 0.75 & 0.1667 & 0.4375 \\ |
| 179 |
> |
H & 0.5735 & 0.2202 & 0.4836 \\ |
| 180 |
> |
H & 0.7420 & 0.0517 & 0.4836 \\ |
| 181 |
> |
O & 0.25 & 0.6667 & 0.4375 \\ |
| 182 |
> |
H & 0.2580 & 0.6693 & 0.3071 \\ |
| 183 |
> |
H & 0.4265 & 0.7255 & 0.4756 \\ \hline |
| 184 |
|
\end{tabular} |
| 185 |
|
\end{wraptable} |
| 186 |
|
|
| 119 |
– |
To construct the basal and prismatic interfaces, the crystallographic |
| 120 |
– |
coordinates above were used to construct an orthorhombic unit cell |
| 121 |
– |
which was then replicated in all three dimensions yielding a |
| 122 |
– |
proton-ordered block of ice I$_{h}$. To expose the desired face, the |
| 123 |
– |
orthorhombic representation was then cut along the ($001$) or ($100$) |
| 124 |
– |
planes for the basal and prismatic faces respectively. The resulting |
| 125 |
– |
block was rotated so that the exposed faces were aligned with the $z$-dimension normal to the exposed face. The block was then cut along two |
| 126 |
– |
perpendicular directions in a way that allowed for perfect periodic |
| 127 |
– |
replication in the $x$ and $y$ axes, creating a slab with either the |
| 128 |
– |
basal or prismatic faces exposed along the $z$ axis. The slab was |
| 129 |
– |
then replicated in the $x$ and $y$ dimensions until a desired sample |
| 130 |
– |
size was obtained. |
| 131 |
– |
|
| 187 |
|
Although experimental solid/liquid coexistant temperature under normal |
| 188 |
|
pressure are close to 273K, Haymet \emph{et al.} have done extensive |
| 189 |
|
work on characterizing the ice/water |
| 190 |
|
interface.\cite{Karim88,Karim90,Hayward01,Bryk02,Hayward02} They have |
| 191 |
< |
found for the SPC/E water model,\cite{Berendsen87} which is also used |
| 192 |
< |
in this study, the ice/water interface is most stable at |
| 193 |
< |
225$\pm$5K.\cite{Bryk02} To create a ice / water interface, a box of |
| 194 |
< |
liquid water that had the same dimensions in $x$ and $y$ was |
| 195 |
< |
equilibrated at 225 K and 1 atm of pressure in the NPAT ensemble (with |
| 196 |
< |
the $z$ axis allowed to fluctuate to equilibrate to the correct |
| 197 |
< |
pressure). The liquid and solid systems were combined by carving out |
| 198 |
< |
any water molecule from the liquid simulation cell that was within 3 |
| 199 |
< |
\AA\ of any atom in the ice slab. |
| 191 |
> |
found that for the SPC/E water model,\cite{Berendsen87} which is also |
| 192 |
> |
used in this study, the ice/water interface is most stable at |
| 193 |
> |
225$\pm$5K.\cite{Bryk02} Therefore, we created our ice / water |
| 194 |
> |
interfaces, utilizing a box of liquid water that had the same |
| 195 |
> |
dimensions in $x$ and $y$ was equilibrated at 225 K and 1 atm of |
| 196 |
> |
pressure in the NPAT ensemble (with the $z$ axis allowed to fluctuate |
| 197 |
> |
to equilibrate to the correct pressure). The liquid and solid systems |
| 198 |
> |
were combined by carving out any water molecule from the liquid |
| 199 |
> |
simulation cell that was within 3 \AA\ of any atom in the ice slab. |
| 200 |
|
|
| 201 |
|
Molecular translation and orientational restraints were applied in the |
| 202 |
|
early stages of equilibration to prevent melting of the ice slab. |
| 203 |
|
These restraints were removed during NVT equilibration, well before |
| 204 |
< |
data collection was carried out. |
| 204 |
> |
data collection was carried out. |
| 205 |
|
|
| 206 |
|
\subsection{Shearing ice / water interfaces without bulk melting} |
| 207 |
|
|
| 208 |
< |
As one drags a solid through a liquid, there will be frictional |
| 209 |
< |
heating that will act to melt the interface. To study the frictional |
| 210 |
< |
behavior of the interface without causing the interface to melt, it is |
| 211 |
< |
necessary to apply a weak thermal gradient along with the momentum |
| 212 |
< |
gradient. This can be accomplished with of the newly-developed |
| 213 |
< |
approaches to reverse non-equilibrium molecular dynamics (RNEMD). The |
| 214 |
< |
velocity shearing and scaling (VSS) variant of RNEMD utilizes a series |
| 215 |
< |
of simultaneous velocity exchanges between two regions within the |
| 216 |
< |
simulation cell.\cite{Kuang12} One of these regions is centered within |
| 217 |
< |
the ice slab, while the other is centrally located in the liquid phase |
| 208 |
> |
As a solid is dragged through a liquid, there is frictional heating |
| 209 |
> |
that will act to melt the interface. To study the behavior of the |
| 210 |
> |
interface under a shear stress without causing the interface to melt, |
| 211 |
> |
it is necessary to apply a weak thermal gradient along with the |
| 212 |
> |
momentum gradient. This can be accomplished using he velocity |
| 213 |
> |
shearing and scaling (VSS) variant of reverse non-equilibrium |
| 214 |
> |
molecular dynamics (RNEMD), which utilizes a series of simultaneous |
| 215 |
> |
velocity exchanges between two regions within the simulation |
| 216 |
> |
cell.\cite{Kuang12} One of these regions is centered within the ice |
| 217 |
> |
slab, while the other is centrally located in the liquid phase |
| 218 |
|
region. VSS-RNEMD provides a set of conservation constraints for |
| 219 |
|
simultaneously creating either a momentum flux or a thermal flux (or |
| 220 |
|
both) between the two slabs. Satisfying the constraint equations |
| 221 |
|
ensures that the new configurations are sampled from the same NVE |
| 222 |
< |
ensemble as previously. |
| 222 |
> |
ensemble as before the VSS move. |
| 223 |
|
|
| 224 |
|
The VSS moves are applied periodically to scale and shift the particle |
| 225 |
|
velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and |
| 274 |
|
|
| 275 |
|
As the simulation progresses, the VSS moves are performed on a regular |
| 276 |
|
basis, and the system develops a thermal and/or velocity gradient in |
| 277 |
< |
response to the applied flux. In a bulk material it is quite simple |
| 277 |
> |
response to the applied flux. In a bulk material, it is quite simple |
| 278 |
|
to use the slope of the temperature or velocity gradients to obtain |
| 279 |
< |
the thermal conductivity or shear viscosity. |
| 279 |
> |
either the thermal conductivity or shear viscosity. |
| 280 |
|
|
| 281 |
|
The VSS-RNEMD approach is versatile in that it may be used to |
| 282 |
|
implement thermal and shear transport simultaneously. Perturbations |
| 286 |
|
map out the shear viscosity of SPC/E water over a wide range of |
| 287 |
|
temperatures (90~K) with a single 1 ns simulation.\cite{Kuang12} |
| 288 |
|
|
| 289 |
< |
Here we are using this method primarily to generate a shear between |
| 290 |
< |
the ice slab and the liquid phase, while using a weak thermal gradient |
| 291 |
< |
to maintaining the interface at the 225K target value. This will |
| 292 |
< |
insure minimal melting of the bulk ice phase and allows us to control |
| 293 |
< |
the exact temperature of the interface. |
| 289 |
> |
For this work, we are using the VSS-RNEMD method primarily to generate |
| 290 |
> |
a shear between the ice slab and the liquid phase, while using a weak |
| 291 |
> |
thermal gradient to maintaining the interface at the 225K target |
| 292 |
> |
value. This will insure minimal melting of the bulk ice phase and |
| 293 |
> |
allows us to control the exact temperature of the interface. |
| 294 |
|
|
| 295 |
|
\subsection{Computational Details} |
| 296 |
< |
All simulations were performed using OpenMD with a time step of 2 fs, |
| 297 |
< |
and periodic boundary conditions in all three dimensions. The systems |
| 298 |
< |
were divided into 100 artificial bins along the $z$-axis for the |
| 299 |
< |
VSS-RNEMD moves, which were attempted every 50 fs. The gradients were |
| 300 |
< |
allowed to develop for 1 ns before data collection was began. Once |
| 301 |
< |
established, four successive 0.5 ns runs were performed for each shear rate. During these simulations, snapshots of the system were taken every 1 ps, and the |
| 247 |
< |
average velocities and densities of each bin were accumulated every |
| 248 |
< |
attempted VSS-RNEMD move. |
| 296 |
> |
All simulations were performed using OpenMD,\cite{OOPSE,openmd} with a |
| 297 |
> |
time step of 2 fs and periodic boundary conditions in all three |
| 298 |
> |
dimensions. Electrostatics were handled using the damped-shifted |
| 299 |
> |
force real-space electrostatic kernel.\cite{Ewald} The systems were |
| 300 |
> |
divided into 100 bins along the $z$-axis for the VSS-RNEMD moves, |
| 301 |
> |
which were attempted every 50 fs. |
| 302 |
|
|
| 303 |
< |
%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. |
| 304 |
< |
%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. |
| 305 |
< |
%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. |
| 306 |
< |
%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. |
| 303 |
> |
The interfaces were equilibrated for a total of 10 ns at equilibrium |
| 304 |
> |
conditions before being exposed to either a shear or thermal gradient. |
| 305 |
> |
This consisted of 5 ns under a constant temperature (NVT) integrator |
| 306 |
> |
set to 225K followed by 5 ns under a microcanonical integrator. Weak |
| 307 |
> |
thermal gradients were allowed to develop using the VSS-RNEMD (NVE) |
| 308 |
> |
integrator using a a small thermal flux ($-2.0\times 10^{-6}$ |
| 309 |
> |
kcal/mol/\AA$^2$/fs) for a duration of 5 ns to allow the gradient to |
| 310 |
> |
stabilize. The resulting temperature gradient was less than 5K over |
| 311 |
> |
the entire 1 nm box length, which was sufficient to keep the |
| 312 |
> |
temperature at the interface within $\pm 1$ K of the 225K target. |
| 313 |
|
|
| 314 |
+ |
Velocity gradients were then imposed using the VSS-RNEMD (NVE) |
| 315 |
+ |
integrator with a range of momentum fluxes. These gradients were |
| 316 |
+ |
allowed to stabilize for 1 ns before data collection began. Once |
| 317 |
+ |
established, four successive 0.5 ns runs were performed for each shear |
| 318 |
+ |
rate. During these simulations, snapshots of the system were taken |
| 319 |
+ |
every 1 ps, and statistics on the structure and dynamics in each bin |
| 320 |
+ |
were accumulated throughout the simulations. |
| 321 |
+ |
|
| 322 |
|
\section{Results and discussion} |
| 323 |
|
|
| 324 |
|
\subsection{Measuring the Width of the Interface} |
| 325 |
< |
\subsubsection{Tetrahedrality Order Parameter} |
| 326 |
< |
Any parameter or function that varies across the interface from a bulk liquid value to a solid value can be used as a measure of the width of the interface. However, due to the VSS-RNEMD moves pertrurbing the momentum of the molecules, parameters such as the translational order parameter and the diffusion order parameter may be artifically skewed. A structural parameter such as the pairwise correlation function would not be influenced by the perturbations. Here, the local order tetraherdal parameter as described by Kumar\cite{Kumar09} and |
| 327 |
< |
Errington\cite{Errington01} was used as a measure of the interfacial width. |
| 325 |
> |
Any order parameter or correlation function that varies across the |
| 326 |
> |
interface from a bulk liquid to a solid can be used as a measure of |
| 327 |
> |
the width of the interface. However, because VSS-RNEMD imposes a |
| 328 |
> |
lateral flow, parameters that depend on translational motion of the |
| 329 |
> |
molecules (e.g. the diffusion constant) may be artifically skewed by |
| 330 |
> |
the RNEMD moves. A structural parameter like a radial distribution |
| 331 |
> |
function is not influenced by the RNEMD perturbations to the same |
| 332 |
> |
degree. Here, we have used the local tetraherdal order parameter as |
| 333 |
> |
described by Kumar\cite{Kumar09} and Errington\cite{Errington01} as a |
| 334 |
> |
measure of the interfacial width. |
| 335 |
|
|
| 336 |
< |
The local tetrahedral order parameter, $q$, is given by |
| 336 |
> |
The local tetrahedral order parameter, $q(z)$, is given by |
| 337 |
|
\begin{equation} |
| 338 |
< |
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 |
| 338 |
> |
q(z) \equiv \int_0^L \sum_{k=1}^{N} \Bigg(1 -\frac{3}{8}\sum_{i=1}^{3} |
| 339 |
> |
\sum_{j=i+1}^{4} \bigg[\cos\psi_{ikj}+\frac{1}{3}\bigg]^2\Bigg) |
| 340 |
> |
\delta(z_{k}-z)\mathrm{d}z |
| 341 |
> |
\label{eq:qz} |
| 342 |
|
\end{equation} |
| 343 |
< |
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. |
| 343 |
> |
where $\psi_{ikj}$ is the angle formed between the oxygen site on |
| 344 |
> |
molecule $k$, and the oxygen sites on its two closest neighbors, |
| 345 |
> |
molecules $i$ and $j$. The local tetrahedral order parameter function |
| 346 |
> |
has a range of (0,1), where the larger the value $q$ has the more |
| 347 |
> |
tetrahedral the ordering of the local environment is. A $q$ value of |
| 348 |
> |
one describes a perfectly tetrahedral environment relative to it and |
| 349 |
> |
its four nearest neighbors, and the parameter's value decreases as the |
| 350 |
> |
local ordering becomes less tetrahedral. Equation \ref{eq:qz} |
| 351 |
> |
describes a $z$-binned tetrahedral order parameter in which the $z$ |
| 352 |
> |
coordinate of the central molecule is used to give a spatial |
| 353 |
> |
description of the local orientational ordering. |
| 354 |
|
|
| 355 |
< |
%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. |
| 355 |
> |
The system was divided into 100 bins along the $z$-dimension, and a |
| 356 |
> |
cutoff radius for the neighboring molecules was set to 3.41 \AA\ . |
| 357 |
> |
The $q_{z}$ values for each snapshot were then averaged to give a |
| 358 |
> |
tetrahedrality profile of the system about the $z$-dimension. The |
| 359 |
> |
profile was then fit with a hyperbolic tangent function given by |
| 360 |
|
|
| 270 |
– |
The system was divided into 100 bins of length $L$ along the $z$-dimension, and a cutoff radius for the neighboring molecules was set to 3.41 \AA\ . A $q_{z}$ value was then determined by averaging the $q$ values for each molecule in the bin. |
| 271 |
– |
\begin{equation} |
| 272 |
– |
q_{z} \equiv \int_0^L \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 |
| 273 |
– |
\end{equation} |
| 274 |
– |
The $q_{z}$ values for each snapshot were then averaged to give a tetrahedrality profile of the system about the $z$-dimension. The profile was then fit with a hyperbolic tangent function given by |
| 275 |
– |
|
| 361 |
|
\begin{equation}\label{tet_fit} |
| 362 |
|
q_{z} \approx q_{liq}+\frac{q_{ice}-q_{liq}}{2}\Bigg[\tanh\bigg(\frac{z-I_{L,m}}{w}\bigg)-\tanh\bigg(\frac{z-I_{R,m}}{w}\bigg)\Bigg]+\beta|(z-z_{mid})| |
| 363 |
|
\end{equation} |
| 434 |
|
|
| 435 |
|
|
| 436 |
|
\section{Conclusion} |
| 437 |
< |
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 apparent independence of the shear rate on the interfacial width, which was found to be 3.2$\pm$0.4 \AA\ and 3.6$\pm$0.2 \AA\ for the basal and prismatic systems. The orientational time correlation function was calculated from varying displacements from the interface. Decomposition by a triexponential decay also showed an apparent independence of the shear rate. The short time decay due to the restoring forces of existing hydrogen bonds decreased at close displacements from the interface, while the middle and long time decays were found to increase. There is also an apparent displacement, $d_{basal}$ and $d_{prismatic}$, from the interface at which these deviations from bulk liquid values occurs. We found $d_{basal}$ and $d_{prismatic}$ to be approximately 15 \AA\ and 8 \AA\ . This implies that the dynamics of water molecules which are structurally equivalent to bulk phase molecules are being perturbed by the presence of the ice and/or the interface. The coefficient of friction of each of the facets was also determined. They were found to be (0.07$\pm$0.01) amu fs\textsuperscript{-1} and (0.032$\pm$0.007) amu fs\textsuperscript{-1} for the basal and prismatic facets respectively. We believe the large difference between the two friction coefficients is due to the direction of the shear and the surface structure of the crystal facets. |
| 437 |
> |
Here we have simulated the basal and prismatic facets of an SPC/E model of the ice I$_\mathrm{h}$ / water interface. Using VSS-RNEMD, the ice was sheared relative to the liquid while imposed thermal gradients kept the interface at a stable temperature. Caculation of the local tetrahedrality order parameter has shown an apparent independence of the shear rate on the interfacial width, which was found to be 3.2$\pm$0.4 \AA\ and 3.6$\pm$0.2 \AA\ for the basal and prismatic systems. The orientational time correlation function was calculated from varying displacements from the interface. Decomposition by a triexponential decay also showed an apparent independence of the shear rate. The short time decay due to the restoring forces of existing hydrogen bonds decreased at close displacements from the interface, while the middle and long time decays were found to increase. There is also an apparent displacement, $d_{basal}$ and $d_{prismatic}$, from the interface at which these deviations from bulk liquid values occurs. We found $d_{basal}$ and $d_{prismatic}$ to be approximately 15 \AA\ and 8 \AA\ . This implies that the dynamics of water molecules which are structurally equivalent to bulk phase molecules are being perturbed by the presence of the ice and/or the interface. The coefficient of friction of each of the facets was also determined. They were found to be (0.07$\pm$0.01) amu fs\textsuperscript{-1} and (0.032$\pm$0.007) amu fs\textsuperscript{-1} for the basal and prismatic facets respectively. We believe the large difference between the two friction coefficients is due to the direction of the shear and the surface structure of the crystal facets. |
| 438 |
|
|
| 439 |
|
\begin{acknowledgement} |
| 440 |
|
Support for this project was provided by the National Science |
