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{MullerPlathe97,MullerPlathe99} 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 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).

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{Viscardy2007a,Viscardy2007b}) 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 hot 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 cold bin 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 P_h^\alpha + \alpha^\prime P_c^\alpha
\end{equation}
where
\begin{equation}
\begin{array}{rcl}
P_h^\alpha & = & \sum_{i = 1}^{N_h} m_i \left[\vec{v}_i\right]_\alpha \\
P_c^\alpha & = & \sum_{j = 1}^{N_c} m_j \left[\vec{v}_j\right]_\alpha
\\
\end{array}
\label{eq:momentumdef}
\end{equation}
Therefore, for each of the three directions, the cold scaling
parameters are a simple function of the hot scaling parameters and
the instantaneous linear momentum in each of the two slabs.
\begin{equation}
\alpha^\prime = 1 + (1 - \alpha) \frac{P_h^\alpha}{P_c^\alpha}.
\label{eq:hotcoldscaling}
\end{equation}

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}
\alpha^2 K_h^\alpha + \left(\alpha^\prime\right)^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
cold scaling parameters ($\alpha^\prime$) from
Eq. (\ref{eq:hotcoldscaling}), we obtain the {\it constraint ellipsoid}:
\begin{equation}
C_x (1+x)^2 + C_y (1+y)^2 + C_z (1+z)^2 = 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}
C_\alpha = \left( K_h^\alpha + K_c^\alpha
  \left(\frac{P_h^\alpha}{P_c^\alpha}\right)^2\right) 
\label{eq:constraintEllipsoidConsts}
\end{equation}
This ellipsoid defines the set of hot slab scaling parameters which can be
applied while preserving both linear momentum in all three directions
as well as total 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 increase in kinetic energy of the hot bin:
\begin{equation}
(x^2-1) K_h^x  + (y^2-1) K_h^y + (z^2-1) K_h^z = J_z
\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_h^x + y^2 K_h^y + z^2 K_h^z = (J_z + K_h^x + K_h^y + K_h^z)
\label{eq:fluxEllipsoid}
\end{equation}
The spatial extent of the {\it flux ellipsoid} is governed both by a
targetted value, $J_z$ as well as the instantaneous values of the
kinetic energy components in the hot 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-space.


One may also define momentum flux (say along the x-direction) as:
\begin{equation}
(x-1) P_h^x  = j_z(p_x)
\end{equation}


\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:	3524
Committed:	Mon Sep 21 20:28:37 2009 UTC (14 years, 9 months ago) by gezelter
Content type:	application/x-tex
File size:	9172 byte(s)
Log Message:	initial import
#	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			the simulation cell.\cite{MullerPlathe97,MullerPlathe99} The
60			artificial flux is typically created by periodically "swapping" either
61			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			directional ($J_z(p_x)$) or isotropic ($J_z$), and the response of a
65			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			J_z(p_x) & = & -\eta \frac{\partial v_x}{\partial z}\\
71			J & = & \lambda \frac{\partial T}{\partial z}
72			\end{eqnarray}
73			RNEMD been widely used to provide computational estimates of thermal
74			conductivities and shear viscosities in a wide range of materials,
75			from liquid copper to monatomic liquids to molecular fluids
76			(e.g. ionic liquids).
77
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			al.}\cite{Viscardy2007a,Viscardy2007b}) because these rely on
86			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			moves. For molecules $\{i\}$ located within the hot slab,
125			\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			located in the cold bin will see a concomitant scaling of velocities,
135			\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			P_h^\alpha + P_c^\alpha = \alpha P_h^\alpha + \alpha^\prime P_c^\alpha
148			\end{equation}
149			where
150			\begin{equation}
151			\begin{array}{rcl}
152			P_h^\alpha & = & \sum_{i = 1}^{N_h} m_i \left[\vec{v}_i\right]_\alpha \\
153			P_c^\alpha & = & \sum_{j = 1}^{N_c} m_j \left[\vec{v}_j\right]_\alpha
154			\\
155			\end{array}
156			\label{eq:momentumdef}
157			\end{equation}
158			Therefore, for each of the three directions, the cold scaling
159			parameters are a simple function of the hot scaling parameters and
160			the instantaneous linear momentum in each of the two slabs.
161			\begin{equation}
162			\alpha^\prime = 1 + (1 - \alpha) \frac{P_h^\alpha}{P_c^\alpha}.
163			\label{eq:hotcoldscaling}
164			\end{equation}
165
166			Conservation of total energy also places constraints on the scaling:
167			\begin{equation}
168			\sum_{\alpha = x,y,z} K_h^\alpha + K_c^\alpha = \sum_{\alpha = x,y,z}
169			\alpha^2 K_h^\alpha + \left(\alpha^\prime\right)^2 K_c^\alpha.
170			\end{equation}
171			where the kinetic energies, $K_h^\alpha$ and $K_c^\alpha$, are computed
172			for each of the three directions in a similar manner to the linear momenta
173			(Eq. \ref{eq:momentumdef}). Substituting in the expressions for the
174			cold scaling parameters ($\alpha^\prime$) from
175			Eq. (\ref{eq:hotcoldscaling}), we obtain the {\it constraint ellipsoid}:
176			\begin{equation}
177			C_x (1+x)^2 + C_y (1+y)^2 + C_z (1+z)^2 = 0,
178			\label{eq:constraintEllipsoid}
179			\end{equation}
180			where the constants are obtained from the instantaneous values of the
181			linear momenta and kinetic energies for the hot and cold slabs,
182			\begin{equation}
183			C_\alpha = \left( K_h^\alpha + K_c^\alpha
184			\left(\frac{P_h^\alpha}{P_c^\alpha}\right)^2\right)
185			\label{eq:constraintEllipsoidConsts}
186			\end{equation}
187			This ellipsoid defines the set of hot slab scaling parameters which can be
188			applied while preserving both linear momentum in all three directions
189			as well as total energy.
190
191			The goal of using velocity scaling variables is to transfer linear
192			momentum or kinetic energy from the cold slab to the hot slab. If the
193			hot and cold slabs are separated along the z-axis, the energy flux is
194			given simply by the increase in kinetic energy of the hot bin:
195			\begin{equation}
196			(x^2-1) K_h^x + (y^2-1) K_h^y + (z^2-1) K_h^z = J_z
197			\end{equation}
198			The expression for the energy flux can be re-written as another
199			ellipsoid centered on $(x,y,z) = 0$:
200			\begin{equation}
201			x^2 K_h^x + y^2 K_h^y + z^2 K_h^z = (J_z + K_h^x + K_h^y + K_h^z)
202			\label{eq:fluxEllipsoid}
203			\end{equation}
204			The spatial extent of the {\it flux ellipsoid} is governed both by a
205			targetted value, $J_z$ as well as the instantaneous values of the
206			kinetic energy components in the hot bin.
207
208			To satisfy an energetic flux as well as the conservation constraints,
209			it is sufficient to determine the points ${x,y,z}$ which lie on both
210			the constraint ellipsoid (Eq. \ref{eq:constraintEllipsoid}) and the
211			flux ellipsoid (Eq. \ref{eq:fluxEllipsoid}), i.e. the intersection of
212			the two ellipsoids in 3-space.
213
214
215			One may also define momentum flux (say along the x-direction) as:
216			\begin{equation}
217			(x-1) P_h^x = j_z(p_x)
218			\end{equation}
219
220
221
222
223
224			\section{Acknowledgments}
225			Support for this project was provided by the National Science
226			Foundation under grant CHE-0848243. Computational time was provided by
227			the Center for Research Computing (CRC) at the University of Notre
228			Dame. \newpage
229
230			\bibliographystyle{jcp2}
231			\bibliography{nivsRnemd}
232			\end{doublespace}
233			\end{document}
234