--- 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)\\