trunk/iceWater/iceWater.tex

\documentclass[%
 jcp,
 jmp,%
 amsmath,amssymb,
 preprint,%
% reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines


\begin{document}

%\preprint{AIP/123-QED}

\title[VSS-RNEMD Simulation of Shearing the Basal and Prismatic Faces of Ice I_{h}]{VSS-RNEMD Simulation of Shearing the Basal and Prismatic Faces of Ice I_{h}}% Force line breaks with \\
%\thanks{Footnote to title of article.}

\author{P. B. Louden}
 %\altaffiliation[Also at]{Physics Department, XYZ University.}%Lines break automatically or can be forced with \\
\author{J. D. Gezelter}
 \homepage{To whom correspondence should be addressed. E-mail: gezelter@nd.edu}
\affiliation{%
Department of Chemistry and Biochemistry, The University of Notre Dame, Notre Dame, IN 46556%\\This line break forced% with \\
}%

\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\maketitle

\section{Abstract}
We have investigated the structural properties of the basal and prismatic facets of an SPC/E model of the ice Ih / water interface when the solid phase is being drawn through liquid water (i.e. sheared relative to the fluid phase). To impose the shear, we utilized a reverse non-equilibrium molecular dynamics (RNEMD) method that creates non-equilibrium conditions using velocity shearing and scaling (VSS) moves of the molecules in two physically separated slabs in the simulation cell. This method can create simultaneous temperature and velocity gradients and allow the measurement of friction and thermal transport properties at interfaces. We present calculations of the interfacial friction coefficients and the apparent independence of shear rate on interfacial width and show that water moving over a flat ice/water interface is close to the no-slip boundary condition.


\section{Introduction}

%Other people looking at the ice/water interface 
%Geologists are concerned with the flow of water over ice
%Antifreeze protein in fish--Haymet's group has cited this before


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.


\section{Methodology}
\subsection{System Construction}
To construct the basal and prismatic systems, first the ice lattices were created. Hirsch and Ojam\"{a}e recently determined possible proton-ordered structures of ice Ih for an orthorhombic unit cell containing eight water molecules. \cite{Hirsch04} The crystallographic coordinates for structure 6 (P$2_{1}2_{1}2_{1}$) were used to construct an orthorhombic unit cell which was then replicated in all three dimensions yielding a proton-ordered block of ice Ih. To expose the desired face, the ice block was then cut along the ($0001$) plane or the ($10\overline{1}0$) plane for the basal and prismatic faces respectively. The ice block was also cut perpendicular to the initial cuts, and oriented so that the desired face is exposed to the $z$-axis. The ice block was then replicated in the $x$ and $y$ dimensions, and lastly liquid phase water molecules were added to the system. Haymet \emph{et al.} have done extensive work on studying and characterizing the ice/water interface. They have found for the SPC/E water model\cite{Berendsen87} (used here), the ice/water interface is most stable at 225$\pm$5K. Therefore, the average temperature of each simulation was 225K. Molecular translation and orientation resrtaints were imposed in the early stages of equilibration to prevent melting of the ice block. These restraints were removed during NVT equilibration, long before data collection. 

\subsection{Computational Details}
All simulations were performed using OpenMD with a time step of 2 fs, and periodic boundary conditions in all three dimensions. The systems were divided into 100 artificial bins along the $z$-axis for the VSS-RNEMD moves, which were attempted every 50 fs. The gradients were allowed to develop for 1 ns before data collection was began. Once established,  snapshots of the system were taken every 1 ps, and the average velocities and densities of each bin were accumulated every attempted VSS-RNEMD move.


%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. 
%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.
%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.
%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.

\subsection{Measuring the Width of the Interface}
In order to characterize the ice/water interface, the local tetrahedral order parameter as described by Kumar\cite{Kumar09} and Errinton\cite{Errington01} was used. The local tetrahedral order parameter, $q$, is given by 
\begin{equation}
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
\end{equation}
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 function decreases as the local ordering becomes less tetrahedral.
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.
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 
\begin{equation}
q_{z}=q_{liq}+\frac{q_{ice}-q_{liq}}{2}(\tanh(\alpha(z-I_{L,m}))-\tanh(\alpha(z-I_{R,m})))
\end{equation}
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,
\begin{equation}
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})|
\end{equation}
where $\beta$ is a fitting parameter and $z_{mid}$ is the midpoint of the z dimension of the simulation box. 


\section{Results and discussion}

%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.

\subsection{Interfacial Width}
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. The interfacial width as described by the fit of the tetrahedrally profile shows no dependence on shear rate as seen in Figure \ref{fig:FIGURENAME}.

\subsection{Coefficient of Friction of the Interface}
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{Balasumramnian99}
\begin{equation}
%<F_{x}^{w}>_{NE}(t)=-S\lambda_{wall}v_{x}(y_{wall})
\langle F_{x}^{w}\rangle(t)=-S\lambda_{wall}v_{x}(y_{wall})
\end{equation}
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 
\begin{equation}
J_{z}(p_{x})=-\lambda_{ice}v_{x}(y_{ice}).
\end{equation}


\section{Conclusion}
Write your conclusion here.

\section{Acknowledgements}
I would like to acknowledge my roommate for not strangling me, yet.

\section{References}
\bibliography{iceWater}%may need to be in same folder?

%Below is an example of how a Figure is imported and referenced in the paper.
%here, FIGURENAME = fig:dip
%\pagebreak
%\begin{figure}
%\includegraphics{dipolemoment.eps}% Here is how to import EPS art
%\caption{\label{fig:dip} The condensation coefficient as a function of the dipole moment.}
%\end{figure}


\end{document}
Revision:	3898
Committed:	Wed Jun 19 17:35:30 2013 UTC (11 years ago) by plouden
Content type:	application/x-tex
File size:	10319 byte(s)
Log Message:	changes
#	User	Rev	Content
1	gezelter	3897	\documentclass[%
2			jcp,
3			jmp,%
4			amsmath,amssymb,
5			preprint,%
6			% reprint,%
7			%author-year,%
8			%author-numerical,%
9			]{revtex4-1}
10
11			\usepackage{graphicx}% Include figure files
12			\usepackage{dcolumn}% Align table columns on decimal point
13			\usepackage{bm}% bold math
14			%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
15			%\linenumbers\relax % Commence numbering lines
16
17
18			\begin{document}
19
20			%\preprint{AIP/123-QED}
21
22			\title[VSS-RNEMD Simulation of Shearing the Basal and Prismatic Faces of Ice I_{h}]{VSS-RNEMD Simulation of Shearing the Basal and Prismatic Faces of Ice I_{h}}% Force line breaks with \\
23			%\thanks{Footnote to title of article.}
24
25			\author{P. B. Louden}
26			%\altaffiliation[Also at]{Physics Department, XYZ University.}%Lines break automatically or can be forced with \\
27			\author{J. D. Gezelter}
28			\homepage{To whom correspondence should be addressed. E-mail: gezelter@nd.edu}
29			\affiliation{%
30			Department of Chemistry and Biochemistry, The University of Notre Dame, Notre Dame, IN 46556%\\This line break forced% with \\
31			}%
32
33			\date{\today}% It is always \today, today,
34			% but any date may be explicitly specified
35
36			\maketitle
37
38			\section{Abstract}
39			We have investigated the structural properties of the basal and prismatic facets of an SPC/E model of the ice Ih / water interface when the solid phase is being drawn through liquid water (i.e. sheared relative to the fluid phase). To impose the shear, we utilized a reverse non-equilibrium molecular dynamics (RNEMD) method that creates non-equilibrium conditions using velocity shearing and scaling (VSS) moves of the molecules in two physically separated slabs in the simulation cell. This method can create simultaneous temperature and velocity gradients and allow the measurement of friction and thermal transport properties at interfaces. We present calculations of the interfacial friction coefficients and the apparent independence of shear rate on interfacial width and show that water moving over a flat ice/water interface is close to the no-slip boundary condition.
40
41
42			\section{Introduction}
43
44			%Other people looking at the ice/water interface
45			%Geologists are concerned with the flow of water over ice
46			%Antifreeze protein in fish--Haymet's group has cited this before
47
48	plouden	3898
49
50	gezelter	3897	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.
51
52
53	plouden	3898
54	gezelter	3897	\section{Methodology}
55			\subsection{System Construction}
56	plouden	3898	To construct the basal and prismatic systems, first the ice lattices were created. Hirsch and Ojam\"{a}e recently determined possible proton-ordered structures of ice Ih for an orthorhombic unit cell containing eight water molecules. \cite{Hirsch04} The crystallographic coordinates for structure 6 (P$2_{1}2_{1}2_{1}$) were used to construct an orthorhombic unit cell which was then replicated in all three dimensions yielding a proton-ordered block of ice Ih. To expose the desired face, the ice block was then cut along the ($0001$) plane or the ($10\overline{1}0$) plane for the basal and prismatic faces respectively. The ice block was also cut perpendicular to the initial cuts, and oriented so that the desired face is exposed to the $z$-axis. The ice block was then replicated in the $x$ and $y$ dimensions, and lastly liquid phase water molecules were added to the system. Haymet \emph{et al.} have done extensive work on studying and characterizing the ice/water interface. They have found for the SPC/E water model\cite{Berendsen87} (used here), the ice/water interface is most stable at 225$\pm$5K. Therefore, the average temperature of each simulation was 225K. Molecular translation and orientation resrtaints were imposed in the early stages of equilibration to prevent melting of the ice block. These restraints were removed during NVT equilibration, long before data collection.
57	gezelter	3897
58			\subsection{Computational Details}
59			All simulations were performed using OpenMD with a time step of 2 fs, and periodic boundary conditions in all three dimensions. The systems were divided into 100 artificial bins along the $z$-axis for the VSS-RNEMD moves, which were attempted every 50 fs. The gradients were allowed to develop for 1 ns before data collection was began. Once established, snapshots of the system were taken every 1 ps, and the average velocities and densities of each bin were accumulated every attempted VSS-RNEMD move.
60
61
62			%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.
63			%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.
64			%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.
65			%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.
66
67			\subsection{Measuring the Width of the Interface}
68			In order to characterize the ice/water interface, the local tetrahedral order parameter as described by Kumar\cite{Kumar09} and Errinton\cite{Errington01} was used. The local tetrahedral order parameter, $q$, is given by
69			\begin{equation}
70			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
71			\end{equation}
72			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 function decreases as the local ordering becomes less tetrahedral.
73			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.
74			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
75			\begin{equation}
76			q_{z}=q_{liq}+\frac{q_{ice}-q_{liq}}{2}(\tanh(\alpha(z-I_{L,m}))-\tanh(\alpha(z-I_{R,m})))
77			\end{equation}
78			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,
79			\begin{equation}
80			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})\|
81			\end{equation}
82			where $\beta$ is a fitting parameter and $z_{mid}$ is the midpoint of the z dimension of the simulation box.
83
84
85			\section{Results and discussion}
86	plouden	3898
87			%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.
88
89	gezelter	3897	\subsection{Interfacial Width}
90			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. The interfacial width as described by the fit of the tetrahedrally profile shows no dependence on shear rate as seen in Figure \ref{fig:FIGURENAME}.
91
92			\subsection{Coefficient of Friction of the Interface}
93			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{Balasumramnian99}
94			\begin{equation}
95			%<F_{x}^{w}>_{NE}(t)=-S\lambda_{wall}v_{x}(y_{wall})
96			\langle F_{x}^{w}\rangle(t)=-S\lambda_{wall}v_{x}(y_{wall})
97			\end{equation}
98			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
99			\begin{equation}
100			J_{z}(p_{x})=-\lambda_{ice}v_{x}(y_{ice}).
101			\end{equation}
102
103
104			\section{Conclusion}
105			Write your conclusion here.
106
107			\section{Acknowledgements}
108			I would like to acknowledge my roommate for not strangling me, yet.
109
110			\section{References}
111			\bibliography{iceWater}%may need to be in same folder?
112
113			%Below is an example of how a Figure is imported and referenced in the paper.
114			%here, FIGURENAME = fig:dip
115			%\pagebreak
116			%\begin{figure}
117			%\includegraphics{dipolemoment.eps}% Here is how to import EPS art
118			%\caption{\label{fig:dip} The condensation coefficient as a function of the dipole moment.}
119			%\end{figure}
120
121
122			\end{document}