trunk/nivsRnemd/nivsRnemd.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\usepackage{tabularx}
\usepackage{graphicx}
%\usepackage{booktabs}
%\usepackage{bibentry}
%\usepackage{mathrsfs}
\usepackage[ref]{overcite}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

\renewcommand\citemid{\ } % no comma in optional referenc note

\begin{document}

\title{A gentler approach to RNEMD: Non-isotropic Velocity Scaling for computing Thermal Conductivity and Shear Viscosity}

\author{Shenyu Kuang and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Introduction}
The original formulation of Reverse Non-equilibrium Molecular Dynamics
(RNEMD) obtains transport coefficients (thermal conductivity and shear
viscosity) in a fluid by imposing an artificial momentum flux between
two thin parallel slabs of material that are spatially separated in
the simulation cell.\cite{MullerPlathe:1997xw,ISI:000080382700030} The
artificial flux is typically created by periodically ``swapping'' either
the entire momentum vector $\vec{p}$ or single components of this
vector ($p_x$) between molecules in each of the two slabs.  If the two
slabs are separated along the z coordinate, the imposed flux is either
directional ($j_z(p_x)$) or isotropic ($J_z$), and the response of a
simulated system to the imposed momentum flux will typically be a
velocity or thermal gradient.  The transport coefficients (shear
viscosity and thermal conductivity) are easily obtained by assuming
linear response of the system,
\begin{eqnarray}
j_z(p_x) & = & -\eta \frac{\partial v_x}{\partial z}\\
J & = & \lambda \frac{\partial T}{\partial z}
\end{eqnarray}
RNEMD has been widely used to provide computational estimates of thermal
conductivities and shear viscosities in a wide range of materials,
from liquid copper to monatomic liquids to molecular fluids
(e.g. ionic liquids).\cite{ISI:000246190100032}

RNEMD is preferable in many ways to the forward NEMD methods because
it imposes what is typically difficult to measure (a flux or stress)
and it is typically much easier to compute momentum gradients or
strains (the response).  For similar reasons, RNEMD is also preferable
to slowly-converging equilibrium methods for measuring thermal
conductivity and shear viscosity (using Green-Kubo relations or the
Helfand moment approach of Viscardy {\it et
  al.}\cite{Viscardy:2007bh,Viscardy:2007lq}) because these rely on
computing difficult to measure quantities.

Another attractive feature of RNEMD is that it conserves both total
linear momentum and total energy during the swaps (as long as the two
molecules have the same identity), so the swapped configurations are
typically samples from the same manifold of states in the
microcanonical ensemble.

Recently, Tenney and Maginn have discovered some problems with the
original RNEMD swap technique.  Notably, large momentum fluxes
(equivalent to frequent momentum swaps between the slabs) can result
in "notched", "peaked" and generally non-thermal momentum
distributions in the two slabs, as well as non-linear thermal and
velocity distributions along the direction of the imposed flux ($z$).
Tenney and Maginn obtained reasonable limits on imposed flux and
self-adjusting metrics for retaining the usability of the method.

In this paper, we develop and test a method for non-isotropic velocity
scaling (NIVS-RNEMD) which retains the desirable features of RNEMD
(conservation of linear momentum and total energy, compatibility with
periodic boundary conditions) while establishing true thermal
distributions in each of the two slabs.  In the next section, we
develop the method for determining the scaling constraints.  We then
test the method on both single component, multi-component, and
non-isotropic mixtures and show that it is capable of providing
reasonable estimates of the thermal conductivity and shear viscosity
in these cases.

\section{Methodology}
We retain the basic idea of Muller-Plathe's RNEMD method; the periodic
system is partitioned into a series of thin slabs along a particular
axis ($z$).  One of the slabs at the end of the periodic box is
designated the ``hot'' slab, while the slab in the center of the box
is designated the ``cold'' slab.  The artificial momentum flux will be
established by transferring momentum from the cold slab and into the
hot slab.

Rather than using momentum swaps, we use a series of velocity scaling
moves.  For molecules $\{i\}$ located within the cold slab,
\begin{equation}
\vec{v}_i \leftarrow \left( \begin{array}{c}
x \\
y \\
z \\
\end{array} \right) \cdot \vec{v}_i
\end{equation}
where ${x, y, z}$ are a set of 3 scaling variables for each of the
three directions in the system.  Likewise, the molecules $\{j\}$
located in the hot slab will see a concomitant scaling of velocities,
\begin{equation}
\vec{v}_j \leftarrow \left( \begin{array}{c}
x^\prime \\
y^\prime \\
z^\prime \\
\end{array} \right) \cdot \vec{v}_j
\end{equation}

Conservation of linear momentum in each of the three directions
($\alpha = x,y,z$) ties the values of the hot and cold bin scaling
parameters together:
\begin{equation}
P_h^\alpha + P_c^\alpha = \alpha^\prime P_h^\alpha + \alpha P_c^\alpha
\end{equation}
where
\begin{equation}
\begin{array}{rcl}
P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i \left[\vec{v}_i\right]_\alpha \\
P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j \left[\vec{v}_j\right]_\alpha \\
\end{array}
\label{eq:momentumdef}
\end{equation}
Therefore, for each of the three directions, the hot scaling
parameters are a simple function of the cold scaling parameters and
the instantaneous linear momentum in each of the two slabs.
\begin{equation}
\alpha^\prime = 1 + (1 - \alpha) p_\alpha
\label{eq:hotcoldscaling}
\end{equation}
where
\begin{equation}
p_\alpha = \frac{P_c^\alpha}{P_h^\alpha}
\end{equation}
for convenience.

Conservation of total energy also places constraints on the scaling:
\begin{equation}
\sum_{\alpha = x,y,z} K_h^\alpha + K_c^\alpha = \sum_{\alpha = x,y,z}
\left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha.
\end{equation}
where the kinetic energies, $K_h^\alpha$ and $K_c^\alpha$, are computed
for each of the three directions in a similar manner to the linear momenta
(Eq. \ref{eq:momentumdef}).  Substituting in the expressions for the
hot scaling parameters ($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}),
we obtain the {\it constraint ellipsoid equation}:
\begin{equation}
\sum_{\alpha = x,y,z} a_\alpha \alpha^2 + b_\alpha \alpha + c_\alpha = 0,
\label{eq:constraintEllipsoid}
\end{equation}
where the constants are obtained from the instantaneous values of the
linear momenta and kinetic energies for the hot and cold slabs,
\begin{equation}
\begin{array}{rcl}
a_\alpha & = &\left(K_c^\alpha + K_h^\alpha
  \left(p_\alpha\right)^2\right) \\
b_\alpha & = & -2 K_h^\alpha p_\alpha \left( 1 + p_\alpha\right) \\
c_\alpha & = & K_h^\alpha p_\alpha^2 + 2 K_h^\alpha p_\alpha - K_c^\alpha \\
\end{array}
\label{eq:constraintEllipsoidConsts}
\end{equation}
This ellipsoid equation defines the set of cold slab scaling
parameters which can be applied while preserving both linear momentum
in all three directions as well as kinetic energy.

The goal of using velocity scaling variables is to transfer linear
momentum or kinetic energy from the cold slab to the hot slab.  If the
hot and cold slabs are separated along the z-axis, the energy flux is
given simply by the decrease in kinetic energy of the cold bin:
\begin{equation}
(1-x^2) K_c^x  + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z\Delta t
\end{equation}
The expression for the energy flux can be re-written as another
ellipsoid centered on $(x,y,z) = 0$:
\begin{equation}
x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = K_c^x + K_c^y + K_c^z - J_z\Delta t
\label{eq:fluxEllipsoid}
\end{equation}
The spatial extent of the {\it flux ellipsoid equation} is governed
both by a targetted value, $J_z$ as well as the instantaneous values of the
kinetic energy components in the cold bin.

To satisfy an energetic flux as well as the conservation constraints,
it is sufficient to determine the points ${x,y,z}$ which lie on both
the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the
flux ellipsoid (Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of
the two ellipsoids in 3-dimensional space.

One may also define momentum flux (say along the x-direction) as:
\begin{equation}
(1-x) P_c^x  = j_z(p_x)\Delta t
\label{eq:fluxPlane}
\end{equation}
The above {\it flux equation} is essentially a plane which is
perpendicular to the x-axis, with its position governed both by a
targetted value, $j_z(p_x)$ as well as the instantaneous value of the
momentum along the x-direction.

Similarly, to satisfy a momentum flux as well as the conservation
constraints, it is sufficient to determine the points ${x,y,z}$ which
lie on both the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid})
and the flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of
an ellipsoid and a plane in 3-dimensional space.

To summarize, by solving respective equation sets, one can determine
possible sets of scaling variables for cold slab. And corresponding
sets of scaling variables for hot slab can be determine as well.

The following problem will be choosing an optimal set of scaling
variables among the possible sets. Although this method is inherently
non-isotropic, the goal is still to maintain the system as isotropic
as possible. Under this consideration, one would like the kinetic
energies in different directions could become as close as each other
after each scaling. Simultaneously, one would also like each scaling
as gentle as possible, i.e. ${x,y,z \rightarrow 1}$, in order to avoid
large perturbation to the system. Therefore, one approach to obtain the
scaling variables would be constructing an criteria function, with
constraints as above equation sets, and solving the function's minimum
by method like Lagrange multipliers.

In order to save computation time, we have a different approach to a
relatively good set of scaling variables with much less calculation
than above. Here is the detail of our simplification of the problem.

In the case of kinetic energy transfer, we impose another constraint
${x = y}$, into the equation sets. Consequently, there are two
variables left. And now one only needs to solve a set of two {\it
  ellipses equations}. This problem would be transformed into solving
one quartic equation for one of the two variables. There are known
generic methods that solve real roots of quartic equations. Then one
can determine the other variable and obtain sets of scaling
variables. Among these sets, one can apply the above criteria to
choose the best set, while much faster with only a few sets to choose.

In the case of momentum flux transfer, we impose another constraint to
set the kinetic energy transfer as zero. In another word, we apply
Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$. After that, with one
variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar set
of equations on the above kinetic energy transfer problem. Therefore,
an approach similar to the above would be sufficient for this as well.

\section{Computational Details}
Our simulation consists of a series of systems. All of these
simulations were run with the OOPSE simulation software
package\cite{Meineke:2005gd} integrated with RNEMD methods.

A Lennard-Jones fluid system was built and tested first. In order to
compare our method with swapping RNEMD, a series of simulations were
performed to calculate the shear viscosity and thermal conductivity of
argon. 2592 atoms were in a orthorhombic cell, which was ${10.06\sigma
  \times 10.06\sigma \times 30.18\sigma}$ by size. The reduced density
${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled direct
comparison between our results and others. These simulations used
Verlet time-stepping algorithm with reduced timestep ${\tau^* =
  4.6\times10^{-4}}$. 

For shear viscosity calculation, the reduced temperature was ${T^* =
  k_B T/\varepsilon = 0.72}$. Simulations were run in microcanonical
ensemble (NVE). For the swapping part, Muller-Plathe's algorithm was
adopted.\cite{ISI:000080382700030} The simulation box was under
periodic boundary condition, and devided into ${N = 20}$ slabs. In each swap,
the top slab ${(n = 1)}$ exchange the most negative $x$ momentum with the
most positive $x$ momentum in the center slab ${(n = N/2 + 1)}$. Referred
to Tenney {\it et al.}\cite{tenneyANDmaginn}, a series of swapping
frequency were chosen. According to each result from swapping
RNEMD, scaling RNEMD simulations were run with the target momentum
flux set to produce a similar momentum flux and shear
rate. Furthermore, various scaling frequencies can be tested for one
single swapping rate. To compare the performance between swapping and
scaling method, temperatures of different dimensions in all the slabs
were observed. 

After each simulation, the shear viscosity was calculated in reduced
unit. The momentum flux was calculated with total unphysical
transferred momentum ${P_x}$ and simulation time $t$:
\begin{equation}
j_z(p_x) = \frac{P_x}{2 t L_x L_y}
\end{equation}
And the velocity gradient ${\langle \partial v_x /\partial z \rangle}$
can be obtained by a linear regression of the velocity profile. From
the shear viscosity $\eta$ calculated with the above parameters, one
can further convert it into reduced unit ${\eta^* = \eta \sigma^2
  (\varepsilon  m)^{-1/2}}$.

For thermal conductivity calculation, simulations were first run under
reduced temperature ${T^* = 0.72}$ in NVE ensemble. Muller-Plathe's
algorithm was adopted in the swapping method. Under identical
simulation box, in each swap, the top slab exchange the molecule with
least kinetic energy with the molecule in the center slab with most
kinetic energy, unless this ``coldest'' molecule in the ``hot'' slab
is still ``hotter'' than the ``hottest'' molecule in the ``cold''
slab. According to swapping RNEMD results, target energy flux for
scaling RNEMD simulations can be set. Also, various scaling
frequencies can be tested for one target energy flux. To compare the
performance between swapping and scaling method, distributions of
velocity and speed in different slabs were observed.

For each swapping rate, thermal conductivity was calculated in reduced
unit. The energy flux was calculated similarly to the momentum flux,
with total unphysical transferred energy ${E_{total}}$ and simulation
time $t$:
\begin{equation}
J_z = \frac{E_{total}}{2 t L_x L_y}
\end{equation}
And the temperature gradient ${\langle\partial T/\partial z\rangle}$
can be obtained by a linear regression of the temperature
profile. From the thermal conductivity $\lambda$ calculated, one can
further convert it into reduced unit ${\lambda^*=\lambda \sigma^2
  m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$.

\section{Results And Discussion}
\subsection{Shear Viscosity}

\section{Acknowledgments}
Support for this project was provided by the National Science
Foundation under grant CHE-0848243. Computational time was provided by
the Center for Research Computing (CRC) at the University of Notre
Dame.  \newpage

\bibliographystyle{jcp2}
\bibliography{nivsRnemd}
\end{doublespace}
\end{document}

Revision:	3534
Committed:	Thu Oct 8 21:42:53 2009 UTC (14 years, 8 months ago) by skuang
Content type:	application/x-tex
File size:	16025 byte(s)
Log Message:	continue simulation details
#	User	Rev	Content
1	gezelter	3524	\documentclass[11pt]{article}
2			\usepackage{amsmath}
3			\usepackage{amssymb}
4			\usepackage{setspace}
5			\usepackage{endfloat}
6			\usepackage{caption}
7			%\usepackage{tabularx}
8			\usepackage{graphicx}
9			%\usepackage{booktabs}
10			%\usepackage{bibentry}
11			%\usepackage{mathrsfs}
12			\usepackage[ref]{overcite}
13			\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
14			\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
15			9.0in \textwidth 6.5in \brokenpenalty=10000
16
17			% double space list of tables and figures
18			\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}
19			\setlength{\abovecaptionskip}{20 pt}
20			\setlength{\belowcaptionskip}{30 pt}
21
22			\renewcommand\citemid{\ } % no comma in optional referenc note
23
24			\begin{document}
25
26			\title{A gentler approach to RNEMD: Non-isotropic Velocity Scaling for computing Thermal Conductivity and Shear Viscosity}
27
28			\author{Shenyu Kuang and J. Daniel
29			Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
30			Department of Chemistry and Biochemistry,\\
31			University of Notre Dame\\
32			Notre Dame, Indiana 46556}
33
34			\date{\today}
35
36			\maketitle
37
38			\begin{doublespace}
39
40			\begin{abstract}
41
42			\end{abstract}
43
44			\newpage
45
46			%\narrowtext
47
48			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49			% BODY OF TEXT
50			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
51
52
53
54			\section{Introduction}
55			The original formulation of Reverse Non-equilibrium Molecular Dynamics
56			(RNEMD) obtains transport coefficients (thermal conductivity and shear
57			viscosity) in a fluid by imposing an artificial momentum flux between
58			two thin parallel slabs of material that are spatially separated in
59	skuang	3534	the simulation cell.\cite{MullerPlathe:1997xw,ISI:000080382700030} The
60	skuang	3531	artificial flux is typically created by periodically ``swapping'' either
61	gezelter	3524	the entire momentum vector $\vec{p}$ or single components of this
62			vector ($p_x$) between molecules in each of the two slabs. If the two
63			slabs are separated along the z coordinate, the imposed flux is either
64	skuang	3532	directional ($j_z(p_x)$) or isotropic ($J_z$), and the response of a
65	gezelter	3524	simulated system to the imposed momentum flux will typically be a
66			velocity or thermal gradient. The transport coefficients (shear
67			viscosity and thermal conductivity) are easily obtained by assuming
68			linear response of the system,
69			\begin{eqnarray}
70	skuang	3532	j_z(p_x) & = & -\eta \frac{\partial v_x}{\partial z}\\
71	gezelter	3524	J & = & \lambda \frac{\partial T}{\partial z}
72			\end{eqnarray}
73	skuang	3528	RNEMD has been widely used to provide computational estimates of thermal
74	gezelter	3524	conductivities and shear viscosities in a wide range of materials,
75			from liquid copper to monatomic liquids to molecular fluids
76	skuang	3528	(e.g. ionic liquids).\cite{ISI:000246190100032}
77	gezelter	3524
78			RNEMD is preferable in many ways to the forward NEMD methods because
79			it imposes what is typically difficult to measure (a flux or stress)
80			and it is typically much easier to compute momentum gradients or
81			strains (the response). For similar reasons, RNEMD is also preferable
82			to slowly-converging equilibrium methods for measuring thermal
83			conductivity and shear viscosity (using Green-Kubo relations or the
84			Helfand moment approach of Viscardy {\it et
85	skuang	3527	al.}\cite{Viscardy:2007bh,Viscardy:2007lq}) because these rely on
86	gezelter	3524	computing difficult to measure quantities.
87
88			Another attractive feature of RNEMD is that it conserves both total
89			linear momentum and total energy during the swaps (as long as the two
90			molecules have the same identity), so the swapped configurations are
91			typically samples from the same manifold of states in the
92			microcanonical ensemble.
93
94			Recently, Tenney and Maginn have discovered some problems with the
95			original RNEMD swap technique. Notably, large momentum fluxes
96			(equivalent to frequent momentum swaps between the slabs) can result
97			in "notched", "peaked" and generally non-thermal momentum
98			distributions in the two slabs, as well as non-linear thermal and
99			velocity distributions along the direction of the imposed flux ($z$).
100			Tenney and Maginn obtained reasonable limits on imposed flux and
101			self-adjusting metrics for retaining the usability of the method.
102
103			In this paper, we develop and test a method for non-isotropic velocity
104			scaling (NIVS-RNEMD) which retains the desirable features of RNEMD
105			(conservation of linear momentum and total energy, compatibility with
106			periodic boundary conditions) while establishing true thermal
107			distributions in each of the two slabs. In the next section, we
108			develop the method for determining the scaling constraints. We then
109			test the method on both single component, multi-component, and
110			non-isotropic mixtures and show that it is capable of providing
111			reasonable estimates of the thermal conductivity and shear viscosity
112			in these cases.
113
114			\section{Methodology}
115			We retain the basic idea of Muller-Plathe's RNEMD method; the periodic
116			system is partitioned into a series of thin slabs along a particular
117			axis ($z$). One of the slabs at the end of the periodic box is
118			designated the ``hot'' slab, while the slab in the center of the box
119			is designated the ``cold'' slab. The artificial momentum flux will be
120			established by transferring momentum from the cold slab and into the
121			hot slab.
122
123			Rather than using momentum swaps, we use a series of velocity scaling
124	skuang	3528	moves. For molecules $\{i\}$ located within the cold slab,
125	gezelter	3524	\begin{equation}
126			\vec{v}_i \leftarrow \left( \begin{array}{c}
127			x \\
128			y \\
129			z \\
130			\end{array} \right) \cdot \vec{v}_i
131			\end{equation}
132			where ${x, y, z}$ are a set of 3 scaling variables for each of the
133			three directions in the system. Likewise, the molecules $\{j\}$
134	skuang	3528	located in the hot slab will see a concomitant scaling of velocities,
135	gezelter	3524	\begin{equation}
136			\vec{v}_j \leftarrow \left( \begin{array}{c}
137			x^\prime \\
138			y^\prime \\
139			z^\prime \\
140			\end{array} \right) \cdot \vec{v}_j
141			\end{equation}
142
143			Conservation of linear momentum in each of the three directions
144			($\alpha = x,y,z$) ties the values of the hot and cold bin scaling
145			parameters together:
146			\begin{equation}
147	skuang	3528	P_h^\alpha + P_c^\alpha = \alpha^\prime P_h^\alpha + \alpha P_c^\alpha
148	gezelter	3524	\end{equation}
149			where
150			\begin{equation}
151			\begin{array}{rcl}
152	skuang	3528	P_c^\alpha & = & \sum_{i = 1}^{N_c} m_i \left[\vec{v}_i\right]_\alpha \\
153			P_h^\alpha & = & \sum_{j = 1}^{N_h} m_j \left[\vec{v}_j\right]_\alpha \\
154	gezelter	3524	\end{array}
155			\label{eq:momentumdef}
156			\end{equation}
157	skuang	3528	Therefore, for each of the three directions, the hot scaling
158			parameters are a simple function of the cold scaling parameters and
159	gezelter	3524	the instantaneous linear momentum in each of the two slabs.
160			\begin{equation}
161	skuang	3528	\alpha^\prime = 1 + (1 - \alpha) p_\alpha
162	gezelter	3524	\label{eq:hotcoldscaling}
163			\end{equation}
164	skuang	3528	where
165			\begin{equation}
166			p_\alpha = \frac{P_c^\alpha}{P_h^\alpha}
167			\end{equation}
168			for convenience.
169	gezelter	3524
170			Conservation of total energy also places constraints on the scaling:
171			\begin{equation}
172			\sum_{\alpha = x,y,z} K_h^\alpha + K_c^\alpha = \sum_{\alpha = x,y,z}
173	skuang	3528	\left(\alpha^\prime\right)^2 K_h^\alpha + \alpha^2 K_c^\alpha.
174	gezelter	3524	\end{equation}
175			where the kinetic energies, $K_h^\alpha$ and $K_c^\alpha$, are computed
176			for each of the three directions in a similar manner to the linear momenta
177			(Eq. \ref{eq:momentumdef}). Substituting in the expressions for the
178	skuang	3528	hot scaling parameters ($\alpha^\prime$) from Eq. (\ref{eq:hotcoldscaling}),
179			we obtain the {\it constraint ellipsoid equation}:
180	gezelter	3524	\begin{equation}
181	skuang	3528	\sum_{\alpha = x,y,z} a_\alpha \alpha^2 + b_\alpha \alpha + c_\alpha = 0,
182	gezelter	3524	\label{eq:constraintEllipsoid}
183			\end{equation}
184			where the constants are obtained from the instantaneous values of the
185			linear momenta and kinetic energies for the hot and cold slabs,
186			\begin{equation}
187	skuang	3528	\begin{array}{rcl}
188			a_\alpha & = &\left(K_c^\alpha + K_h^\alpha
189			\left(p_\alpha\right)^2\right) \\
190			b_\alpha & = & -2 K_h^\alpha p_\alpha \left( 1 + p_\alpha\right) \\
191			c_\alpha & = & K_h^\alpha p_\alpha^2 + 2 K_h^\alpha p_\alpha - K_c^\alpha \\
192			\end{array}
193	gezelter	3524	\label{eq:constraintEllipsoidConsts}
194			\end{equation}
195	skuang	3528	This ellipsoid equation defines the set of cold slab scaling
196			parameters which can be applied while preserving both linear momentum
197	skuang	3530	in all three directions as well as kinetic energy.
198	gezelter	3524
199			The goal of using velocity scaling variables is to transfer linear
200			momentum or kinetic energy from the cold slab to the hot slab. If the
201			hot and cold slabs are separated along the z-axis, the energy flux is
202	skuang	3528	given simply by the decrease in kinetic energy of the cold bin:
203	gezelter	3524	\begin{equation}
204	skuang	3534	(1-x^2) K_c^x + (1-y^2) K_c^y + (1-z^2) K_c^z = J_z\Delta t
205	gezelter	3524	\end{equation}
206			The expression for the energy flux can be re-written as another
207			ellipsoid centered on $(x,y,z) = 0$:
208			\begin{equation}
209	skuang	3534	x^2 K_c^x + y^2 K_c^y + z^2 K_c^z = K_c^x + K_c^y + K_c^z - J_z\Delta t
210	gezelter	3524	\label{eq:fluxEllipsoid}
211			\end{equation}
212	skuang	3529	The spatial extent of the {\it flux ellipsoid equation} is governed
213			both by a targetted value, $J_z$ as well as the instantaneous values of the
214	skuang	3530	kinetic energy components in the cold bin.
215	gezelter	3524
216			To satisfy an energetic flux as well as the conservation constraints,
217			it is sufficient to determine the points ${x,y,z}$ which lie on both
218			the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the
219			flux ellipsoid (Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of
220	skuang	3528	the two ellipsoids in 3-dimensional space.
221	gezelter	3524
222			One may also define momentum flux (say along the x-direction) as:
223			\begin{equation}
224	skuang	3534	(1-x) P_c^x = j_z(p_x)\Delta t
225	skuang	3531	\label{eq:fluxPlane}
226	gezelter	3524	\end{equation}
227	skuang	3531	The above {\it flux equation} is essentially a plane which is
228			perpendicular to the x-axis, with its position governed both by a
229			targetted value, $j_z(p_x)$ as well as the instantaneous value of the
230			momentum along the x-direction.
231	gezelter	3524
232	skuang	3531	Similarly, to satisfy a momentum flux as well as the conservation
233			constraints, it is sufficient to determine the points ${x,y,z}$ which
234			lie on both the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid})
235			and the flux plane (Eq. \ref{eq:fluxPlane}), i.e. the intersection of
236			an ellipsoid and a plane in 3-dimensional space.
237	gezelter	3524
238	skuang	3531	To summarize, by solving respective equation sets, one can determine
239			possible sets of scaling variables for cold slab. And corresponding
240			sets of scaling variables for hot slab can be determine as well.
241	gezelter	3524
242	skuang	3531	The following problem will be choosing an optimal set of scaling
243			variables among the possible sets. Although this method is inherently
244			non-isotropic, the goal is still to maintain the system as isotropic
245			as possible. Under this consideration, one would like the kinetic
246			energies in different directions could become as close as each other
247			after each scaling. Simultaneously, one would also like each scaling
248			as gentle as possible, i.e. ${x,y,z \rightarrow 1}$, in order to avoid
249			large perturbation to the system. Therefore, one approach to obtain the
250			scaling variables would be constructing an criteria function, with
251			constraints as above equation sets, and solving the function's minimum
252			by method like Lagrange multipliers.
253	gezelter	3524
254	skuang	3531	In order to save computation time, we have a different approach to a
255			relatively good set of scaling variables with much less calculation
256			than above. Here is the detail of our simplification of the problem.
257	gezelter	3524
258	skuang	3531	In the case of kinetic energy transfer, we impose another constraint
259			${x = y}$, into the equation sets. Consequently, there are two
260			variables left. And now one only needs to solve a set of two {\it
261			ellipses equations}. This problem would be transformed into solving
262			one quartic equation for one of the two variables. There are known
263			generic methods that solve real roots of quartic equations. Then one
264			can determine the other variable and obtain sets of scaling
265			variables. Among these sets, one can apply the above criteria to
266			choose the best set, while much faster with only a few sets to choose.
267
268			In the case of momentum flux transfer, we impose another constraint to
269			set the kinetic energy transfer as zero. In another word, we apply
270			Eq. \ref{eq:fluxEllipsoid} and let ${J_z = 0}$. After that, with one
271			variable fixed by Eq. \ref{eq:fluxPlane}, one now have a similar set
272			of equations on the above kinetic energy transfer problem. Therefore,
273			an approach similar to the above would be sufficient for this as well.
274
275			\section{Computational Details}
276	skuang	3534	Our simulation consists of a series of systems. All of these
277			simulations were run with the OOPSE simulation software
278			package\cite{Meineke:2005gd} integrated with RNEMD methods.
279	skuang	3531
280	skuang	3532	A Lennard-Jones fluid system was built and tested first. In order to
281			compare our method with swapping RNEMD, a series of simulations were
282			performed to calculate the shear viscosity and thermal conductivity of
283	skuang	3534	argon. 2592 atoms were in a orthorhombic cell, which was ${10.06\sigma
284			\times 10.06\sigma \times 30.18\sigma}$ by size. The reduced density
285			${\rho^* = \rho\sigma^3}$ was thus 0.85, which enabled direct
286			comparison between our results and others. These simulations used
287			Verlet time-stepping algorithm with reduced timestep ${\tau^* =
288			4.6\times10^{-4}}$.
289	skuang	3532
290			For shear viscosity calculation, the reduced temperature was ${T^* =
291	skuang	3534	k_B T/\varepsilon = 0.72}$. Simulations were run in microcanonical
292	skuang	3532	ensemble (NVE). For the swapping part, Muller-Plathe's algorithm was
293			adopted.\cite{ISI:000080382700030} The simulation box was under
294	skuang	3534	periodic boundary condition, and devided into ${N = 20}$ slabs. In each swap,
295			the top slab ${(n = 1)}$ exchange the most negative $x$ momentum with the
296			most positive $x$ momentum in the center slab ${(n = N/2 + 1)}$. Referred
297	skuang	3532	to Tenney {\it et al.}\cite{tenneyANDmaginn}, a series of swapping
298	skuang	3534	frequency were chosen. According to each result from swapping
299	skuang	3532	RNEMD, scaling RNEMD simulations were run with the target momentum
300	skuang	3534	flux set to produce a similar momentum flux and shear
301			rate. Furthermore, various scaling frequencies can be tested for one
302			single swapping rate. To compare the performance between swapping and
303			scaling method, temperatures of different dimensions in all the slabs
304			were observed.
305	skuang	3532
306	skuang	3534	After each simulation, the shear viscosity was calculated in reduced
307			unit. The momentum flux was calculated with total unphysical
308			transferred momentum ${P_x}$ and simulation time $t$:
309			\begin{equation}
310			j_z(p_x) = \frac{P_x}{2 t L_x L_y}
311			\end{equation}
312			And the velocity gradient ${\langle \partial v_x /\partial z \rangle}$
313			can be obtained by a linear regression of the velocity profile. From
314			the shear viscosity $\eta$ calculated with the above parameters, one
315			can further convert it into reduced unit ${\eta^* = \eta \sigma^2
316			(\varepsilon m)^{-1/2}}$.
317	skuang	3532
318	skuang	3534	For thermal conductivity calculation, simulations were first run under
319			reduced temperature ${T^* = 0.72}$ in NVE ensemble. Muller-Plathe's
320			algorithm was adopted in the swapping method. Under identical
321			simulation box, in each swap, the top slab exchange the molecule with
322			least kinetic energy with the molecule in the center slab with most
323			kinetic energy, unless this ``coldest'' molecule in the ``hot'' slab
324			is still ``hotter'' than the ``hottest'' molecule in the ``cold''
325			slab. According to swapping RNEMD results, target energy flux for
326			scaling RNEMD simulations can be set. Also, various scaling
327			frequencies can be tested for one target energy flux. To compare the
328			performance between swapping and scaling method, distributions of
329			velocity and speed in different slabs were observed.
330
331			For each swapping rate, thermal conductivity was calculated in reduced
332			unit. The energy flux was calculated similarly to the momentum flux,
333			with total unphysical transferred energy ${E_{total}}$ and simulation
334			time $t$:
335			\begin{equation}
336			J_z = \frac{E_{total}}{2 t L_x L_y}
337			\end{equation}
338			And the temperature gradient ${\langle\partial T/\partial z\rangle}$
339			can be obtained by a linear regression of the temperature
340			profile. From the thermal conductivity $\lambda$ calculated, one can
341			further convert it into reduced unit ${\lambda^*=\lambda \sigma^2
342			m^{1/2} k_B^{-1}\varepsilon^{-1/2}}$.
343
344			\section{Results And Discussion}
345			\subsection{Shear Viscosity}
346
347	gezelter	3524	\section{Acknowledgments}
348			Support for this project was provided by the National Science
349			Foundation under grant CHE-0848243. Computational time was provided by
350			the Center for Research Computing (CRC) at the University of Notre
351			Dame. \newpage
352
353			\bibliographystyle{jcp2}
354			\bibliography{nivsRnemd}
355			\end{doublespace}
356			\end{document}
357