ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/iceWater/iceWater.tex
Revision: 3922
Committed: Tue Jul 16 19:15:15 2013 UTC (10 years, 11 months ago) by plouden
Content type: application/x-tex
File size: 26419 byte(s)
Log Message:
fixed q_z equations

File Contents

# User Rev Content
1 gezelter 3914 \documentclass[journal = jpccck, manuscript = article]{achemso}
2     \setkeys{acs}{usetitle = true}
3     \usepackage{achemso}
4     \usepackage{natbib}
5     \usepackage{multirow}
6     \usepackage{wrapfig}
7     \usepackage{fixltx2e}
8     %\mciteErrorOnUnknownfalse
9 gezelter 3897
10 gezelter 3914 \usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions
11     \usepackage{url}
12 gezelter 3897
13    
14 gezelter 3914 \title{Do the facets of ice $I_\mathrm{h}$ crystals have different
15     friction coefficients? Simulating shear in ice/water interfaces}
16 gezelter 3897
17     \author{P. B. Louden}
18 gezelter 3914 \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}
23 gezelter 3897
24 gezelter 3914 \keywords{}
25 gezelter 3897
26 gezelter 3914 \begin{document}
27 gezelter 3897
28 gezelter 3914 \begin{abstract}
29     We have investigated the structural properties of the basal and
30     prismatic facets of an SPC/E model of the ice Ih / water interface
31     when the solid phase is being drawn through liquid water (i.e. sheared
32     relative to the fluid phase). To impose the shear, we utilized a
33     reverse non-equilibrium molecular dynamics (RNEMD) method that creates
34     non-equilibrium conditions using velocity shearing and scaling (VSS)
35     moves of the molecules in two physically separated slabs in the
36     simulation cell. This method can create simultaneous temperature and
37     velocity gradients and allow the measurement of friction transport
38     properties at interfaces. We present calculations of the interfacial
39     friction coefficients and the apparent independence of shear rate on
40     interfacial width and show that water moving over a flat ice/water
41     interface is close to the no-slip boundary condition.
42     \end{abstract}
43 gezelter 3897
44 gezelter 3914 \newpage
45 gezelter 3897
46     \section{Introduction}
47 plouden 3919 %-----Outline of Intro---------------
48     % in general, ice/water interface is important b/c ....
49     % here are some people who have worked on ice/water, trying to understand the processes above ....
50     % with the recent development of VSS-RNEMD, we can now look at the shearing problem
51     % talk about what we will present in this paper
52     % -------End Intro------------------
53 gezelter 3897
54 plouden 3919 %Gay02: cites many other ice/water papers, make sure to cite them.
55    
56 plouden 3920 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}.
57 plouden 3919
58 plouden 3920 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.
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. 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.
60 plouden 3919
61 gezelter 3897 \section{Methodology}
62    
63 gezelter 3914 \subsection{Stable ice I$_\mathrm{h}$ / water interfaces}
64    
65     The structure of ice I$_\mathrm{h}$ is well understood; it
66     crystallizes in a hexagonal space group P$6_3/mmc$, and the hexagonal
67     crystals of ice have two faces that are commonly exposed, the basal
68     face $\{0~0~0~1\}$, which forms the top and bottom of each hexagonal
69     plate, and the prismatic $\{1~0~\bar{1}~0\}$ face which forms the
70     sides of the plate. Other less-common, but still important, faces of
71     ice I$_\mathrm{h}$ are the secondary prism face, $\{1~1~\bar{2}~0\}$,
72     and the prismatic face, $\{2~0~\bar{2}~1\}$.
73    
74     Ice I$_\mathrm{h}$ is normally proton disordered in bulk crystals,
75     although the surfaces probably have a preference for proton ordering
76     along strips of dangling H-atoms and Oxygen lone
77     pairs.\cite{Buch:2008fk}
78    
79     For small simulated ice interfaces, it is useful to have a
80     proton-ordered, but zero-dipole crystal that exposes these strips of
81     dangling H-atoms and lone pairs. Also, if we're going to place
82     another material in contact with one of the ice crystalline planes, it
83     is useful to have an orthorhombic (rectangular) box to work with. A
84     recent paper by Hirsch and Ojam\"{a}e describes how to create
85     proton-ordered bulk ice I$_\mathrm{h}$ in alternative orthorhombic
86 gezelter 3915 cells.\cite{Hirsch04}
87 gezelter 3914
88     We have using Hirsch and Ojam\"{a}e's structure 6 which is an
89     orthorhombic cell ($P2_12_12_1$) that produces a proton-ordered
90     version of ice Ih. Table \ref{tab:equiv} contains a mapping between
91     the Miller indices in the P$6_3/mmc$ crystal system and those in the
92     Hirsch and Ojam\"{a}e $P2_12_12_1$ system.
93    
94     \begin{wraptable}{r}{3.5in}
95     \begin{tabular}{|ccc|} \hline
96     & hexagonal & orthorhombic \\
97     & ($P6_3/mmc$) & ($P2_12_12_1$) \\
98     crystal face & Miller indices & equivalent \\ \hline
99     basal & $\{0~0~0~1\}$ & $\{0~0~1\}$ \\
100     prism & $\{1~0~\bar{1}~0\}$ & $\{1~0~0\}$ \\
101     secondary prism & $\{1~1~\bar{2}~0\}$ & $\{1~3~0\}$ \\
102     pryamidal & $\{2~0~\bar{2}~1\}$ & $\{2~0~1\}$ \\ \hline
103     \end{tabular}
104     \end{wraptable}
105    
106     Structure 6 from the Hirsch and Ojam\"{a}e paper has lattice
107     parameters $a = 4.49225$, $b = 7.78080$, $c = 7.33581$ and two water
108     molecules whose atoms reside at the following fractional coordinates:
109    
110     \begin{wraptable}{r}{3.25in}
111     \begin{tabular}{|ccccc|} \hline
112     atom label & type & x & y & z \\ \hline
113     O$_{a}$ & O & 0.75 & 0.1667 & 0.4375 \\
114     H$_{1a}$ & H & 0.5735 & 0.2202 & 0.4836 \\
115     H$_{2a}$ & H & 0.7420 & 0.0517 & 0.4836 \\
116     O$_{b}$ & O & 0.25 & 0.6667 & 0.4375 \\
117     H$_{1b}$ & H & 0.2580 & 0.6693 & 0.3071 \\
118     H$_{2b}$ & H & 0.4265 & 0.7255 & 0.4756 \\ \hline
119     \end{tabular}
120     \end{wraptable}
121    
122     To construct the basal and prismatic interfaces, the crystallographic
123     coordinates above were used to construct an orthorhombic unit cell
124     which was then replicated in all three dimensions yielding a
125     proton-ordered block of ice I$_{h}$. To expose the desired face, the
126     orthorhombic representation was then cut along the ($001$) or ($100$)
127     planes for the basal and prismatic faces respectively. The resulting
128     block was rotated so that the exposed faces were aligned with the $z$
129     axis normal to the exposed face. The block was then cut along two
130     perpendicular directions in a way that allowed for perfect periodic
131     replication in the $x$ and $y$ axes, creating a slab with either the
132     basal or prismatic faces exposed along the $z$ axis. The slab was
133     then replicated in the $x$ and $y$ dimensions until a desired sample
134     size was obtained.
135    
136     Although experimental solid/liquid coexistant temperature under normal
137     pressure are close to 273K, Haymet \emph{et al.} have done extensive
138     work on characterizing the ice/water
139     interface.\cite{Karim88,Karim90,Hayward01,Bryk02,Hayward02} They have
140     found for the SPC/E water model,\cite{Berendsen87} which is also used
141     in this study, the ice/water interface is most stable at
142     225$\pm$5K.\cite{Bryk02} To create a ice / water interface, a box of
143     liquid water that had the same dimensions in $x$ and $y$ was
144     equilibrated at 225 K and 1 atm of pressure in the NPAT ensemble (with
145     the $z$ axis allowed to fluctuate to equilibrate to the correct
146     pressure). The liquid and solid systems were combined by carving out
147     any water molecule from the liquid simulation cell that was within 3
148     \AA\ of any atom in the ice slab.
149    
150     Molecular translation and orientational restraints were applied in the
151     early stages of equilibration to prevent melting of the ice slab.
152     These restraints were removed during NVT equilibration, well before
153     data collection was carried out.
154    
155 gezelter 3918 \subsection{Shearing ice / water interfaces without bulk melting}
156    
157     As one drags a solid through a liquid, there will be frictional
158     heating that will act to melt the interface. To study the frictional
159     behavior of the interface without causing the interface to melt, it is
160     necessary to apply a weak thermal gradient along with the momentum
161     gradient. This can be accomplished with of the newly-developed
162     approaches to reverse non-equilibrium molecular dynamics (RNEMD). The
163     velocity shearing and scaling (VSS) variant of RNEMD utilizes a series
164     of simultaneous velocity exchanges between two regions within the
165     simulation cell.\cite{Kuang12} One of these regions is centered within
166     the ice slab, while the other is centrally located in the liquid phase
167     region. VSS-RNEMD provides a set of conservation constraints for
168     simultaneously creating either a momentum flux or a thermal flux (or
169     both) between the two slabs. Satisfying the constraint equations
170     ensures that the new configurations are sampled from the same NVE
171     ensemble as previously.
172    
173     The VSS moves are applied periodically to scale and shift the particle
174     velocities ($\mathbf{v}_i$ and $\mathbf{v}_j$) in two slabs ($H$ and
175     $C$) which are separated by half of the simulation box,
176     \begin{displaymath}
177     \begin{array}{rclcl}
178    
179     & \underline{\mathrm{shearing}} & &
180     \underline{~~~~~~~~~~~~\mathrm{scaling}~~~~~~~~~~~~} \\
181     \mathbf{v}_i \leftarrow &
182     \mathbf{a}_c\ & + & c\cdot\left(\mathbf{v}_i - \langle\mathbf{v}_c
183     \rangle\right) + \langle\mathbf{v}_c\rangle \\
184     \mathbf{v}_j \leftarrow &
185     \mathbf{a}_h & + & h\cdot\left(\mathbf{v}_j - \langle\mathbf{v}_h
186     \rangle\right) + \langle\mathbf{v}_h\rangle .
187    
188     \end{array}
189     \end{displaymath}
190     Here $\langle\mathbf{v}_c\rangle$ and $\langle\mathbf{v}_h\rangle$ are
191     the center of mass velocities in the $C$ and $H$ slabs, respectively.
192     Within the two slabs, particles receive incremental changes or a
193     ``shear'' to their velocities. The amount of shear is governed by the
194     imposed momentum flux, $\mathbf{j}_z(\mathbf{p})$
195     \begin{eqnarray}
196     \mathbf{a}_c & = & - \mathbf{j}_z(\mathbf{p}) \Delta t / M_c \label{vss1}\\
197     \mathbf{a}_h & = & + \mathbf{j}_z(\mathbf{p}) \Delta t / M_h \label{vss2}
198     \end{eqnarray}
199     where $M_{\{c,h\}}$ is the total mass of particles within each of the
200     slabs and $\Delta t$ is the interval between two separate operations.
201    
202     To simultaneously impose a thermal flux ($J_z$) between the slabs we
203     use energy conservation constraints,
204     \begin{eqnarray}
205     K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\mathbf{v}_c
206     \rangle^2) + \frac{1}{2}M_c (\langle \mathbf{v}_c \rangle + \mathbf{a}_c)^2 \label{vss3}\\
207     K_h + J_z\Delta t & = & h^2 (K_h - \frac{1}{2}M_h \langle\mathbf{v}_h
208     \rangle^2) + \frac{1}{2}M_h (\langle \mathbf{v}_h \rangle +
209     \mathbf{a}_h)^2 \label{vss4}.
210     \label{constraint}
211     \end{eqnarray}
212     Simultaneous solution of these quadratic formulae for the scaling
213     coefficients, $c$ and $h$, will ensure that the simulation samples from
214     the original microcanonical (NVE) ensemble. Here $K_{\{c,h\}}$ is the
215     instantaneous translational kinetic energy of each slab. At each time
216     interval, it is a simple matter to solve for $c$, $h$, $\mathbf{a}_c$,
217     and $\mathbf{a}_h$, subject to the imposed momentum flux,
218     $j_z(\mathbf{p})$, and thermal flux, $J_z$, values. Since the VSS
219     operations do not change the kinetic energy due to orientational
220     degrees of freedom or the potential energy of a system, configurations
221     after the VSS move have exactly the same energy (and linear
222     momentum) as before the move.
223    
224     As the simulation progresses, the VSS moves are performed on a regular
225     basis, and the system develops a thermal and/or velocity gradient in
226     response to the applied flux. In a bulk material it is quite simple
227     to use the slope of the temperature or velocity gradients to obtain
228     the thermal conductivity or shear viscosity.
229    
230     The VSS-RNEMD approach is versatile in that it may be used to
231     implement thermal and shear transport simultaneously. Perturbations
232     of velocities away from the ideal Maxwell-Boltzmann distributions are
233     minimal, as is thermal anisotropy. This ability to generate
234     simultaneous thermal and shear fluxes has been previously utilized to
235     map out the shear viscosity of SPC/E water over a wide range of
236     temperatures (90~K) with a single 1 ns simulation.\cite{Kuang12}
237    
238     Here we are using this method primarily to generate a shear between
239     the ice slab and the liquid phase, while using a weak thermal gradient
240     to maintaining the interface at the 225K target value. This will
241     insure minimal melting of the bulk ice phase and allows us to control
242     the exact temperature of the interface.
243    
244 gezelter 3897 \subsection{Computational Details}
245 gezelter 3918 All simulations were performed using OpenMD with a time step of 2 fs,
246     and periodic boundary conditions in all three dimensions. The systems
247     were divided into 100 artificial bins along the $z$-axis for the
248     VSS-RNEMD moves, which were attempted every 50 fs. The gradients were
249     allowed to develop for 1 ns before data collection was began. Once
250     established, snapshots of the system were taken every 1 ps, and the
251     average velocities and densities of each bin were accumulated every
252     attempted VSS-RNEMD move.
253 gezelter 3897
254     %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.
255     %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.
256     %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.
257     %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.
258    
259     \subsection{Measuring the Width of the Interface}
260 gezelter 3918 In order to characterize the ice/water interface, the local
261     tetrahedral order parameter as described by Kumar\cite{Kumar09} and
262     Errington\cite{Errington01} was used. The local tetrahedral order
263     parameter, $q$, is given by
264 gezelter 3897 \begin{equation}
265     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
266     \end{equation}
267 plouden 3909 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.
268    
269     %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.
270    
271 plouden 3922 The system was divided into 100 bins of length $L$ along the $z$-axis, and a cutoff radius for the neighboring molecules was set to 3.41 \AA\. A $q_{z}$ value was then determined for each bin by averaging the $q$ values for each molecule in the bin.
272 gezelter 3897 \begin{equation}
273 plouden 3921 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
274 gezelter 3897 \end{equation}
275 plouden 3922 The $q_{z}$ values for each snapshot 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
276 plouden 3921 \begin{equation}
277     q_{z} \approx q_{liq}+\frac{q_{ice}-q_{liq}}{2}(\tanh(\alpha(z-I_{L,m}))-\tanh(\alpha(z-I_{R,m})))
278     \end{equation}
279 gezelter 3897 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,
280     \begin{equation}
281 plouden 3921 q_{z} \approx q_{liq}+\frac{q_{ice}-q_{liq}}{2}(\tanh(\alpha(z-I_{L,m}))-\tanh(\alpha(z-I_{R,m})))+\beta|(z-z_{mid})|
282 gezelter 3897 \end{equation}
283     where $\beta$ is a fitting parameter and $z_{mid}$ is the midpoint of the z dimension of the simulation box.
284    
285    
286     \section{Results and discussion}
287 plouden 3898
288     %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.
289    
290 plouden 3909 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.
291 plouden 3904
292 gezelter 3914 \begin{figure}
293     \includegraphics[width=\linewidth]{bComicStrip}
294     \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.}
295     \end{figure}
296 plouden 3904
297 gezelter 3914 %(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.
298    
299     \begin{figure}
300     \includegraphics[width=\linewidth]{pComicStrip}
301     \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.}
302     \end{figure}
303    
304    
305 gezelter 3897 \subsection{Interfacial Width}
306 plouden 3909 %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.
307     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.
308 gezelter 3897
309 plouden 3919 %Need to reword the following paragraph
310 plouden 3921 Beaglehole and Wilson have measured the ice/water interface to have a thickness approximately 10-20 \AA\ for both the basal and prismatic face of ice by ellipticity measurements \cite{Beaglehole93}. While Haymet \emph{et al.} have reported values that agree with these measurements, our results do not. We believe this arises from the different methods used to measure the interfacial width. Haymet and co-workers use the 10-90 widths of the translational, average density, diffusion, and orientational decay times \cite{Hayward01} to measure the interface, whereas we are using the local tetrahedral order parameter.
311 plouden 3919
312 plouden 3921
313     %Beaglehole and Wilson have measured the ice/water interface to have a thickness approximately 10-20 \AA\ for both the basal and prismatic face of ice by ellipticity measurements \cite{Beaglehole93}. Haymet \emph{et al.} agrees with these measurements, our results do not. We are using a parameter from the hyperbolic tangent fit of the local tetrahedrality order parameter to determine the interfacial width, whereas Haymet and co-workers use the 10-90 widths of the translational, average density, diffusion, and orientational decay times \cite{Hayward01}.
314    
315    
316 gezelter 3897 \subsection{Coefficient of Friction of the Interface}
317 gezelter 3915 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}
318 gezelter 3897 \begin{equation}
319     %<F_{x}^{w}>_{NE}(t)=-S\lambda_{wall}v_{x}(y_{wall})
320     \langle F_{x}^{w}\rangle(t)=-S\lambda_{wall}v_{x}(y_{wall})
321     \end{equation}
322     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
323     \begin{equation}
324     J_{z}(p_{x})=-\lambda_{ice}v_{x}(y_{ice}).
325     \end{equation}
326    
327 plouden 3909 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}.
328 gezelter 3897
329 plouden 3904 %Ask dan about truncating versus rounding the values for lambda.
330 plouden 3907 %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}.
331 plouden 3921
332 gezelter 3914 \begin{figure}
333     \includegraphics[width=\linewidth]{CoeffFric}
334     \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 }
335     \end{figure}
336 plouden 3904
337 gezelter 3897 \section{Conclusion}
338 plouden 3909 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.
339 gezelter 3897
340 plouden 3909
341 gezelter 3914 \begin{acknowledgement}
342     Support for this project was provided by the National Science
343     Foundation under grant CHE-0848243. Computational time was provided
344     by the Center for Research Computing (CRC) at the University of
345     Notre Dame.
346     \end{acknowledgement}
347 gezelter 3897
348 gezelter 3914 \newpage
349     \bibstyle{achemso}
350 plouden 3909 \bibliography{iceWater}
351 gezelter 3897
352 gezelter 3914 \begin{tocentry}
353     \begin{wrapfigure}{l}{0.5\textwidth}
354     \begin{center}
355     \includegraphics[width=\linewidth]{SystemImage.png}
356     \end{center}
357     \end{wrapfigure}
358     An image of our system.
359     \end{tocentry}
360 gezelter 3897
361 gezelter 3914 \end{document}
362 gezelter 3897
363 plouden 3904 % basal: slope=11.02808, error in slope = 0.4489844
364     %prismatic: slope = 19.95948, error in slope = 0.5370894