--- stokes/stokes.tex	2011/12/11 03:22:47	3780
+++ stokes/stokes.tex	2011/12/12 23:39:21	3784
@@ -128,17 +128,19 @@ Similar to the NIVS methodology,\cite{kuang:164101} we
 interfaces.
 
 \section{Methodology}
-Similar to the NIVS methodology,\cite{kuang:164101} we consider a
+Similar to the NIVS method,\cite{kuang:164101} we consider a
 periodic system divided into a series of slabs along a certain axis
 (e.g. $z$). The unphysical thermal and/or momentum flux is designated
-from the center slab to one of the end slabs, and thus the center slab
-would have a lower temperature than the end slab (unless the thermal
-flux is negative). Therefore, the center slab is denoted as ``$c$''
-while the end slab as ``$h$''.
+from the center slab to one of the end slabs, and thus the thermal
+flux results in a lower temperature of the center slab than the end
+slab, and the momentum flux results in negative center slab momentum
+with positive end slab momentum (unless these fluxes are set
+negative). Therefore, the center slab is denoted as ``$c$'', while the
+end slab as ``$h$''.
 
-To impose these fluxes, we periodically apply separate operations to
-velocities of particles {$i$} within the center slab and of particles
-{$j$} within the end slab:
+To impose these fluxes, we periodically apply different set of
+operations on velocities of particles {$i$} within the center slab and
+those of particles {$j$} within the end slab:
 \begin{eqnarray}
 \vec{v}_i & \leftarrow & c\cdot\left(\vec{v}_i - \langle\vec{v}_c
   \rangle\right) + \left(\langle\vec{v}_c\rangle + \vec{a}_c\right) \\
@@ -147,7 +149,7 @@ before an operation occurs. When a momentum flux $\vec
 \end{eqnarray}
 where $\langle\vec{v}_c\rangle$ and $\langle\vec{v}_h\rangle$ denotes
 the instantaneous bulk velocity of slabs $c$ and $h$ respectively
-before an operation occurs. When a momentum flux $\vec{j}_z(\vec{p})$
+before an operation is applied. When a momentum flux $\vec{j}_z(\vec{p})$
 presents, these bulk velocities would have a corresponding change
 ($\vec{a}_c$ and $\vec{a}_h$ respectively) according to Newton's
 second law:
@@ -155,17 +157,18 @@ where
 M_c \vec{a}_c & = & -\vec{j}_z(\vec{p}) \Delta t \\
 M_h \vec{a}_h & = &  \vec{j}_z(\vec{p}) \Delta t
 \end{eqnarray}
-where
+where $M$ denotes total mass of particles within a slab:
 \begin{eqnarray}
 M_c & = & \sum_{i = 1}^{N_c} m_i \\
 M_h & = & \sum_{j = 1}^{N_h} m_j
 \end{eqnarray}
-and $\Delta t$ is the interval between two operations.
+and $\Delta t$ is the interval between two separate operations.
 
-The above operations conserve the linear momentum of a periodic
-system. To satisfy total energy conservation as well as to impose a
-thermal flux $J_z$, one would have
-[SUPPORT INFO MIGHT BE NECESSARY TO PUT EXTRA MATH IN]
+The above operations already conserve the linear momentum of a
+periodic system. To further satisfy total energy conservation as well
+as to impose the thermal flux $J_z$, the following equations are
+included as well:
+[MAY PUT EXTRA MATH IN SUPPORT INFO OR APPENDIX]
 \begin{eqnarray}
 K_c - J_z\Delta t & = & c^2 (K_c - \frac{1}{2}M_c \langle\vec{v}_c
 \rangle^2) + \frac{1}{2}M_c (\langle \vec{v}_c \rangle + \vec{a}_c)^2 \\
@@ -173,27 +176,28 @@ $c$ and $h$ respectively before an operation occurs. T
 \rangle^2) + \frac{1}{2}M_h (\langle \vec{v}_h \rangle + \vec{a}_h)^2
 \end{eqnarray}
 where $K_c$ and $K_h$ denotes translational kinetic energy of slabs
-$c$ and $h$ respectively before an operation occurs. These
+$c$ and $h$ respectively before an operation is applied. These
 translational kinetic energy conservation equations are sufficient to
-ensure total energy conservation, as the operations applied do not
-change the potential energy of a system, given that the potential
+ensure total energy conservation, as the operations applied in our
+method do not change the kinetic energy related to other degrees of
+freedom or the potential energy of a system, given that its potential
 energy does not depend on particle velocity.
 
 The above sets of equations are sufficient to determine the velocity
 scaling coefficients ($c$ and $h$) as well as $\vec{a}_c$ and
-$\vec{a}_h$. Note that two roots of $c$ and $h$ exist
-respectively. However, to avoid dramatic perturbations to a system,
-the positive roots (which are closer to 1) are chosen. Figure
-\ref{method} illustrates the implementation of this algorithm in an
-individual step.
+$\vec{a}_h$. Note that there are two roots respectively for $c$ and
+$h$. However, the positive roots (which are closer to 1) are chosen so
+that the perturbations to a system can be reduced to a minimum. Figure
+\ref{method} illustrates the implementation sketch of this algorithm
+in an individual step.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{method}
 \caption{Illustration of the implementation of the algorithm in a
   single step. Starting from an ideal velocity distribution, the
-  transformation is used to apply both thermal and momentum flux from
-  the ``c'' slab to the ``h'' slab. As the figure shows, the thermal
-  distributions preserve after this operation.}
+  transformation is used to apply the effect of both a thermal and a
+  momentum flux from the ``c'' slab to the ``h'' slab. As the figure
+  shows, thermal distributions can preserve after this operation.}
 \label{method}
 \end{figure}
 
@@ -201,80 +205,85 @@ This approach is more computationaly efficient compare
 thermal and/or momentum flux can be applied and the corresponding
 temperature and/or momentum gradients can be established.
 
-This approach is more computationaly efficient compared to the
-previous NIVS method, in that only quadratic equations are involved,
-while the NIVS method needs to solve a quartic equations. Furthermore,
-the method implements isotropic scaling of velocities in respective
-slabs, unlike the NIVS, where an extra criteria function is necessary
-to choose a set of coefficients that performs the most isotropic
-scaling. More importantly, separating the momentum flux imposing from
-velocity scaling avoids the underlying cause that NIVS produced
-thermal anisotropy when applying a momentum flux.
+Compared to the previous NIVS method, this approach is computationally
+more efficient in that only quadratic equations are involved to
+determine a set of scaling coefficients, while the NIVS method needs
+to solve quartic equations. Furthermore, this method implements
+isotropic scaling of velocities in respective slabs, unlike the NIVS,
+where an extra criteria function is necessary to choose a set of
+coefficients that performs a scaling as isotropic as possible. More
+importantly, separating the means of momentum flux imposing from
+velocity scaling avoids the underlying cause to thermal anisotropy in
+NIVS when applying a momentum flux. And later sections will
+demonstrate that this can improve the performance in shear viscosity
+simulations.
 
-The advantages of the approach over the original momentum swapping
-approach lies in its nature to preserve a Gaussian
-distribution. Because the momentum swapping tends to render a
-nonthermal distribution, when the imposed flux is relatively large,
-diffusion of the neighboring slabs could no longer remedy this effect,
-and nonthermal distributions would be observed. Results in later
-section will illustrate this effect.
+The advantages of this approach over the original momentum swapping
+one lies in its nature of preserve a Maxwell-Boltzmann distribution
+(mathimatically a Gaussian function). However, because the momentum
+swapping steps tend to result in a nonthermal distribution, when an
+imposed flux is relatively large, diffusions from the neighboring
+slabs could no longer remedy this effect, problematic distributions
+would be observed. Results in later section will illustrate this
+effect in more detail.
 
 \section{Computational Details}
 The algorithm has been implemented in our MD simulation code,
 OpenMD\cite{Meineke:2005gd,openmd}. We compare the method with
-previous RNEMD methods or equilibrium MD methods in homogeneous fluids
+previous RNEMD methods or equilibrium MD (EMD) methods in homogeneous fluids
 (Lennard-Jones and SPC/E water). And taking advantage of the method,
 we simulate the interfacial friction of different heterogeneous
 interfaces (gold-organic solvent and gold-SPC/E water and ice-liquid
 water).
 
 \subsection{Simulation Protocols}
-The systems to be investigated are set up in a orthorhombic simulation
-cell with periodic boundary conditions in all three dimensions. The
-$z$ axis of these cells were longer and was set as the gradient axis
-of temperature and/or momentum. Thus the cells were divided into $N$
-slabs along this axis, with various $N$ depending on individual
-system. The $x$ and $y$ axis were usually of the same length in
-homogeneous systems or close to each other where interfaces
-presents. In all cases, before introducing a nonequilibrium method to
-establish steady thermal and/or momentum gradients for further
-measurements and calculations, canonical ensemble with a Nos\'e-Hoover
-thermostat\cite{hoover85} and microcanonical ensemble equilibrations
-were used to prepare systems ready for data
-collections. Isobaric-isothermal equilibrations are performed before
-this for SPC/E water systems to reach normal pressure (1 bar), while
-similar equilibrations are used for interfacial systems to relax the
-surface tensions.
+The systems to be investigated are set up in orthorhombic simulation
+cells with periodic boundary conditions in all three dimensions. The
+$z$ axis of these cells were longer and set as the temperature and/or
+momentum gradient axis. And the cells were evenly divided into $N$
+slabs along this axis, with various $N$ depending on individual
+system. The $x$ and $y$ axis were of the same length in homogeneous
+systems or had length scale close to each other where heterogeneous
+interfaces presents. In all cases, before introducing a nonequilibrium
+method to establish steady thermal and/or momentum gradients for
+further measurements and calculations, canonical ensemble with a
+Nos\'e-Hoover thermostat\cite{hoover85} and microcanonical ensemble
+equilibrations were used before data collections. For SPC/E water
+simulations, isobaric-isothermal equilibrations\cite{melchionna93} are
+performed before the above to reach normal pressure (1 bar); for
+interfacial systems, similar equilibrations are used to relax the
+surface tensions of the $xy$ plane.
 
-While homogeneous fluid systems can be set up with random
-configurations, our interfacial systems needs extra steps to ensure
-the interfaces be established properly for computations. The
+While homogeneous fluid systems can be set up with rather random
+configurations, our interfacial systems needs a series of steps to
+ensure the interfaces be established properly for computations. The
 preparation and equilibration of butanethiol covered gold (111)
 surface and further solvation and equilibration process is described
-as in reference \cite{kuang:AuThl}. 
+in details as in reference \cite{kuang:AuThl}.
 
 As for the ice/liquid water interfaces, the basal surface of ice
 lattice was first constructed. Hirsch {\it et
   al.}\cite{doi:10.1021/jp048434u} explored the energetics of ice
-lattices with different proton orders. We refer to their results and
-choose the configuration of the lowest energy after geometry
-optimization as the unit cells of our ice lattices. Although
-experimental solid/liquid coexistant temperature near normal pressure
-is 273K, Bryk and Haymet's simulations of ice/liquid water interfaces
-with different models suggest that for SPC/E, the most stable
-interface is observed at 225$\pm$5K. Therefore, all our ice/liquid
-water simulations were carried out under 225K. To have extra
-protection of the ice lattice during initial equilibration (when the
-randomly generated liquid phase configuration could release large
-amount of energy in relaxation), a constraint method (REF?) was
-adopted until the high energy configuration was relaxed.
-[MAY ADD A FIGURE HERE FOR BASAL PLANE, MAY INCLUDE PRISM IF POSSIBLE]
+lattices with all possible proton order configurations. We refer to
+their results and choose the configuration of the lowest energy after
+geometry optimization as the unit cell for our ice lattices. Although
+experimental solid/liquid coexistant temperature under normal pressure
+should be close to 273K, Bryk and Haymet's simulations of ice/liquid
+water interfaces with different models suggest that for SPC/E, the
+most stable interface is observed at 225$\pm$5K.\cite{bryk:10258}
+Therefore, our ice/liquid water simulations were carried out at
+225K. To have extra protection of the ice lattice during initial
+equilibration (when the randomly generated liquid phase configuration
+could release large amount of energy in relaxation), restraints were
+applied to the ice lattice to avoid inadvertent melting by the heat
+dissipated from the high enery configurations.
+[MAY ADD A SNAPSHOT FOR BASAL PLANE]
 
 \subsection{Force Field Parameters}
 For comparison of our new method with previous work, we retain our
-force field parameters consistent with the results we will compare
-with. The Lennard-Jones fluid used here for argon , and reduced unit
-results are reported for direct comparison purpose.
+force field parameters consistent with previous simulations. Argon is
+the Lennard-Jones fluid used here, and its results are reported in
+reduced unit for direct comparison purpose.
 
 As for our water simulations, SPC/E model is used throughout this work
 for consistency. Previous work for transport properties of SPC/E water
@@ -289,9 +298,9 @@ The small organic molecules included in our simulation
 Spohr potential was adopted\cite{ISI:000167766600035} to depict
 Au-H$_2$O interactions.
 
-The small organic molecules included in our simulations are the Au
-surface capping agent butanethiol and liquid hexane and toluene. The
-United-Atom
+For our gold/organic liquid interfaces, the small organic molecules
+included in our simulations are the Au surface capping agent
+butanethiol and liquid hexane and toluene. The United-Atom
 models\cite{TraPPE-UA.thiols,TraPPE-UA.alkanes,TraPPE-UA.alkylbenzenes}
 for these components were used in this work for better computational
 efficiency, while maintaining good accuracy. We refer readers to our
@@ -319,8 +328,8 @@ two in the denominator is present for the heat transpo
 \end{equation}
 where $L_x$ and $L_y$ denotes the dimensions of the plane in a
 simulation cell perpendicular to the thermal gradient, and a factor of
-two in the denominator is present for the heat transport occurs in
-both $+z$ and $-z$ directions. The temperature gradient
+two in the denominator is necessary for the heat transport occurs in
+both $+z$ and $-z$ directions. The average temperature gradient
 ${\langle\partial T/\partial z\rangle}$ can be obtained by a linear
 regression of the temperature profile, which is recorded during a
 simulation for each slab in a cell. For Lennard-Jones simulations,
@@ -329,18 +338,19 @@ by switching off $J_z$. One can specify the vector
 
 \subsection{Shear viscosities}
 Alternatively, the method can carry out shear viscosity calculations
-by switching off $J_z$. One can specify the vector
-$\vec{j}_z(\vec{p})$ by choosing the three components
+by specify a momentum flux. In our algorithm, one can specify the
+three components of the flux vector $\vec{j}_z(\vec{p})$
 respectively. For shear viscosity simulations, $j_z(p_z)$ is usually
-set to zero. Although for isotropic systems, the direction of
-$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, the ability
-of arbitarily specifying the vector direction in our method provides
-convenience in anisotropic simulations.
+set to zero. For isotropic systems, the direction of
+$\vec{j}_z(\vec{p})$ on the $xy$ plane should not matter, but the
+ability of arbitarily specifying the vector direction in our method
+could provide convenience in anisotropic simulations.
 
-Similar to thermal conductivity computations, linear response of the
-momentum gradient with respect to the shear stress is assumed, and the
-shear viscosity ($\eta$) can be obtained with the imposed momentum
-flux (e.g. in $x$ direction) and the measured gradient:
+Similar to thermal conductivity computations, for a homogeneous
+system, linear response of the momentum gradient with respect to the
+shear stress is assumed, and the shear viscosity ($\eta$) can be
+obtained with the imposed momentum flux (e.g. in $x$ direction) and
+the measured gradient:
 \begin{equation}
 j_z(p_x) = -\eta \frac{\partial v_x}{\partial z}
 \end{equation}
@@ -349,12 +359,14 @@ the data collection time. Also, the velocity gradient
 j_z(p_x) = \frac{P_x}{2 t L_x L_y}
 \end{equation}
 with $P_x$ being the total non-physical momentum transferred within
-the data collection time. Also, the velocity gradient
+the data collection time. Also, the averaged velocity gradient
 ${\langle\partial v_x/\partial z\rangle}$ can be obtained using linear
-regression of the $x$ component of the mean velocity, $\langle
-v_x\rangle$, in each of the bins. For Lennard-Jones simulations, shear
-viscosities are reported in reduced units
+regression of the $x$ component of the mean velocity ($\langle
+v_x\rangle$) in each of the bins. For Lennard-Jones simulations, shear
+viscosities are also reported in reduced units
 (${\eta^*=\eta\sigma^2(\varepsilon m)^{-1/2}}$).
+
+[COMBINE JZKE W JZPX]
 
 \subsection{Interfacial friction and Slip length}
 While the shear stress results in a velocity gradient within bulk
@@ -428,7 +440,8 @@ calculations with various fluxes in reduced units.
         \multicolumn{2}{c}{$\eta^*$} \\
         \hline
         Swap Interval $t^*$ & Equivalent $J_z^*$ or $j_z^*(p_x)$ &
-        NIVS & This work & Swapping & This work \\
+        NIVS\footnote{Reference \cite{kuang:164101}.} & This work &
+        Swapping & This work \\
         \hline
         0.116 & 0.16  & 7.30(0.10) & 6.25(0.23) & 3.57(0.06) & 3.53(0.16)\\
         0.232 & 0.09  & 6.95(0.09) & 8.02(0.56) & 3.64(0.05) & 3.43(0.17)\\
