ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/langevin/langevin.tex

Comparing trunk/langevin/langevin.tex (file contents):
Revision 3310 by gezelter, Fri Jan 11 22:08:57 2008 UTC vs.
Revision 3367 by gezelter, Thu Mar 13 22:16:01 2008 UTC

#	Line 1 | Line 1
1		%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2	<	%\documentclass[aps,prb,preprint]{revtex4}
2	>	%\documentclass[aps,prb,preprint,amsmath,amssymb,endfloats]{revtex4}
3		\documentclass[11pt]{article}
4	<	\usepackage{endfloat}
5	<	\usepackage{amsmath,bm}
4	>	\usepackage{amsmath}
5		\usepackage{amssymb}
6		\usepackage{times}
7		\usepackage{mathptmx}
8		\usepackage{setspace}
9	<	\usepackage{tabularx}
9	>	\usepackage{endfloat}
10	>	%\usepackage{tabularx}
11		\usepackage{graphicx}
12	<	\usepackage{booktabs}
13	<	\usepackage{bibentry}
14	<	\usepackage{mathrsfs}
12	>	%\usepackage{booktabs}
13	>	%\usepackage{bibentry}
14	>	%\usepackage{mathrsfs}
15		\usepackage[ref]{overcite}
16		\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
17		\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
18		9.0in \textwidth 6.5in \brokenpenalty=10000
19	<	\renewcommand{\baselinestretch}{1.2}
19	>
20		\renewcommand\citemid{\ } % no comma in optional referenc note
21
22		\begin{document}
23
24	<	\title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape }
24	>	\title{Langevin dynamics for rigid bodies of arbitrary shape}
25
26	<	\author{Xiuquan Sun, Teng Lin and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
27	<	gezelter@nd.edu} \\
28	<	Department of Chemistry and Biochemistry\\
26	>	\author{Xiuquan Sun, Teng Lin and J. Daniel
27	>	Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
28	>	Department of Chemistry and Biochemistry,\\
29		University of Notre Dame\\
30		Notre Dame, Indiana 46556}
31
32		\date{\today}
33
34	–	\maketitle \doublespacing
34
35	<	\begin{abstract}
35	>	\maketitle
36
37	+
38	+
39	+	\begin{abstract}
40	+	We present an algorithm for carrying out Langevin dynamics simulations
41	+	on complex rigid bodies by incorporating the hydrodynamic resistance
42	+	tensors for arbitrary shapes into an advanced rotational integration
43	+	scheme.  The integrator gives quantitative agreement with both
44	+	analytic and approximate hydrodynamic theories for a number of model
45	+	rigid bodies, and works well at reproducing the solute dynamical
46	+	properties (diffusion constants, and orientational relaxation times)
47	+	obtained from explicitly-solvated simulations.
48		\end{abstract}
49
50		\newpage
#	Line 45 | Line 55 | Notre Dame, Indiana 46556}
55		%                          BODY OF TEXT
56		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
57
58	+	\begin{doublespace}
59	+
60		\section{Introduction}
61
62		%applications of langevin dynamics
63	<	As alternative to Newtonian dynamics, Langevin dynamics, which
64	<	mimics a simple heat bath with stochastic and dissipative forces,
65	<	has been applied in a variety of studies. The stochastic treatment
66	<	of the solvent enables us to carry out substantially longer time
67	<	simulations. Implicit solvent Langevin dynamics simulations of
68	<	met-enkephalin not only outperform explicit solvent simulations for
69	<	computational efficiency, but also agrees very well with explicit
70	<	solvent simulations for dynamical properties.\cite{Shen2002}
59	<	Recently, applying Langevin dynamics with the UNRES model, Liow and
60	<	his coworkers suggest that protein folding pathways can be possibly
61	<	explored within a reasonable amount of time.\cite{Liwo2005} The
62	<	stochastic nature of the Langevin dynamics also enhances the
63	<	sampling of the system and increases the probability of crossing
64	<	energy barriers.\cite{Banerjee2004, Cui2003} Combining Langevin
65	<	dynamics with Kramers's theory, Klimov and Thirumalai identified
66	<	free-energy barriers by studying the viscosity dependence of the
67	<	protein folding rates.\cite{Klimov1997} In order to account for
68	<	solvent induced interactions missing from implicit solvent model,
69	<	Kaya incorporated desolvation free energy barrier into implicit
70	<	coarse-grained solvent model in protein folding/unfolding studies
71	<	and discovered a higher free energy barrier between the native and
72	<	denatured states. Because of its stability against noise, Langevin
73	<	dynamics is very suitable for studying remagnetization processes in
74	<	various systems.\cite{Palacios1998,Berkov2002,Denisov2003} For
75	<	instance, the oscillation power spectrum of nanoparticles from
76	<	Langevin dynamics simulation has the same peak frequencies for
77	<	different wave vectors, which recovers the property of magnetic
78	<	excitations in small finite structures.\cite{Berkov2005a}
63	>	Langevin dynamics, which mimics a heat bath using both stochastic and
64	>	dissipative forces, has been applied in a variety of situations as an
65	>	alternative to molecular dynamics with explicit solvent molecules.
66	>	The stochastic treatment of the solvent allows the use of simulations
67	>	with substantially longer time and length scales. In general, the
68	>	dynamic and structural properties obtained from Langevin simulations
69	>	agree quite well with similar properties obtained from explicit
70	>	solvent simulations.
71
72	<	%review rigid body dynamics
73	<	Rigid bodies are frequently involved in the modeling of different
74	<	areas, from engineering, physics, to chemistry. For example,
75	<	missiles and vehicle are usually modeled by rigid bodies.  The
76	<	movement of the objects in 3D gaming engine or other physics
77	<	simulator is governed by the rigid body dynamics. In molecular
78	<	simulation, rigid body is used to simplify the model in
87	<	protein-protein docking study{\cite{Gray2003}}.
72	>	Recent examples of the usefulness of Langevin simulations include a
73	>	study of met-enkephalin in which Langevin simulations predicted
74	>	dynamical properties that were largely in agreement with explicit
75	>	solvent simulations.\cite{Shen2002} By applying Langevin dynamics with
76	>	the UNRES model, Liwo and his coworkers suggest that protein folding
77	>	pathways can be explored within a reasonable amount of
78	>	time.\cite{Liwo2005}
79
80	<	It is very important to develop stable and efficient methods to
81	<	integrate the equations of motion for orientational degrees of
82	<	freedom. Euler angles are the natural choice to describe the
83	<	rotational degrees of freedom. However, due to $\frac {1}{sin
84	<	\theta}$ singularities, the numerical integration of corresponding
85	<	equations of these motion is very inefficient and inaccurate.
86	<	Although an alternative integrator using multiple sets of Euler
87	<	angles can overcome this difficulty\cite{Barojas1973}, the
88	<	computational penalty and the loss of angular momentum conservation
89	<	still remain. A singularity-free representation utilizing
99	<	quaternions was developed by Evans in 1977.\cite{Evans1977}
100	<	Unfortunately, this approach used a nonseparable Hamiltonian
101	<	resulting from the quaternion representation, which prevented the
102	<	symplectic algorithm from being utilized. Another different approach
103	<	is to apply holonomic constraints to the atoms belonging to the
104	<	rigid body. Each atom moves independently under the normal forces
105	<	deriving from potential energy and constraint forces which are used
106	<	to guarantee the rigidness. However, due to their iterative nature,
107	<	the SHAKE and Rattle algorithms also converge very slowly when the
108	<	number of constraints increases.\cite{Ryckaert1977, Andersen1983}
80	>	The stochastic nature of Langevin dynamics also enhances the sampling
81	>	of the system and increases the probability of crossing energy
82	>	barriers.\cite{Cui2003,Banerjee2004} Combining Langevin dynamics with
83	>	Kramers' theory, Klimov and Thirumalai identified free-energy
84	>	barriers by studying the viscosity dependence of the protein folding
85	>	rates.\cite{Klimov1997} In order to account for solvent induced
86	>	interactions missing from the implicit solvent model, Kaya
87	>	incorporated a desolvation free energy barrier into protein
88	>	folding/unfolding studies and discovered a higher free energy barrier
89	>	between the native and denatured states.\cite{HuseyinKaya07012005}
90
91	<	A break-through in geometric literature suggests that, in order to
92	<	develop a long-term integration scheme, one should preserve the
93	<	symplectic structure of the propagator. By introducing a conjugate
94	<	momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
95	<	equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
96	<	proposed to evolve the Hamiltonian system in a constraint manifold
97	<	by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
98	<	An alternative method using the quaternion representation was
99	<	developed by Omelyan.\cite{Omelyan1998} However, both of these
119	<	methods are iterative and inefficient. In this section, we descibe a
120	<	symplectic Lie-Poisson integrator for rigid bodies developed by
121	<	Dullweber and his coworkers\cite{Dullweber1997} in depth.
91	>	In typical LD simulations, the friction and random ($f_r$) forces on
92	>	individual atoms are taken from Stokes' law,
93	>	\begin{eqnarray}
94	>	m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + f_r(t) \notag \\
95	>	\langle f_r(t) \rangle & = & 0 \\
96	>	\langle f_r(t) f_r(t') \rangle & = & 2 k_B T \xi m \delta(t - t') \notag
97	>	\end{eqnarray}
98	>	where $\xi \approx 6 \pi \eta \rho$.  Here $\eta$ is the viscosity of the
99	>	implicit solvent, and $\rho$ is the hydrodynamic radius of the atom.
100
101	<	%review langevin/browninan dynamics for arbitrarily shaped rigid body
102	<	Combining Langevin or Brownian dynamics with rigid body dynamics,
103	<	one can study slow processes in biomolecular systems. Modeling DNA
104	<	as a chain of rigid beads, which are subject to harmonic potentials
105	<	as well as excluded volume potentials, Mielke and his coworkers
106	<	discovered rapid superhelical stress generations from the stochastic
107	<	simulation of twin supercoiling DNA with response to induced
108	<	torques.\cite{Mielke2004} Membrane fusion is another key biological
109	<	process which controls a variety of physiological functions, such as
110	<	release of neurotransmitters \textit{etc}. A typical fusion event
111	<	happens on the time scale of a millisecond, which is impractical to
112	<	study using atomistic models with newtonian mechanics. With the help
113	<	of coarse-grained rigid body model and stochastic dynamics, the
114	<	fusion pathways were explored by many
115	<	researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the
116	<	difficulty of numerical integration of anisotropic rotation, most of
117	<	the rigid body models are simply modeled using spheres, cylinders,
118	<	ellipsoids or other regular shapes in stochastic simulations. In an
141	<	effort to account for the diffusion anisotropy of arbitrary
142	<	particles, Fernandes and de la Torre improved the original Brownian
143	<	dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
101	>	The use of rigid substructures,\cite{Chun:2000fj}
102	>	coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunX._jp0762020}
103	>	and ellipsoidal representations of protein side
104	>	chains~\cite{Fogolari:1996lr} has made the use of the Stokes-Einstein
105	>	approximation problematic.  A rigid substructure moves as a single
106	>	unit with orientational as well as translational degrees of freedom.
107	>	This requires a more general treatment of the hydrodynamics than the
108	>	spherical approximation provides.  Also, the atoms involved in a rigid
109	>	or coarse-grained structure have solvent-mediated interactions with
110	>	each other, and these interactions are ignored if all atoms are
111	>	treated as separate spherical particles.  The theory of interactions
112	>	{\it between} bodies moving through a fluid has been developed over
113	>	the past century and has been applied to simulations of Brownian
114	>	motion.\cite{FIXMAN:1986lr,Ramachandran1996}
115	>
116	>	In order to account for the diffusion anisotropy of complex shapes,
117	>	Fernandes and Garc\'{i}a de la Torre improved an earlier Brownian
118	>	dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by
119		incorporating a generalized $6\times6$ diffusion tensor and
120	<	introducing a simple rotation evolution scheme consisting of three
121	<	consecutive rotations.\cite{Fernandes2002} Unfortunately, unexpected
122	<	errors and biases are introduced into the system due to the
123	<	arbitrary order of applying the noncommuting rotation
124	<	operators.\cite{Beard2003} Based on the observation the momentum
125	<	relaxation time is much less than the time step, one may ignore the
126	<	inertia in Brownian dynamics. However, the assumption of zero
127	<	average acceleration is not always true for cooperative motion which
128	<	is common in protein motion. An inertial Brownian dynamics (IBD) was
129	<	proposed to address this issue by adding an inertial correction
130	<	term.\cite{Beard2000} As a complement to IBD which has a lower bound
131	<	in time step because of the inertial relaxation time, long-time-step
132	<	inertial dynamics (LTID) can be used to investigate the inertial
133	<	behavior of the polymer segments in low friction
134	<	regime.\cite{Beard2000} LTID can also deal with the rotational
135	<	dynamics for nonskew bodies without translation-rotation coupling by
136	<	separating the translation and rotation motion and taking advantage
137	<	of the analytical solution of hydrodynamics properties. However,
138	<	typical nonskew bodies like cylinders and ellipsoids are inadequate
139	<	to represent most complex macromolecule assemblies. These intricate
140	<	molecules have been represented by a set of beads and their
141	<	hydrodynamic properties can be calculated using variants on the
142	<	standard hydrodynamic interaction tensors.
120	>	introducing a rotational evolution scheme consisting of three
121	>	consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are
122	>	introduced into the system due to the arbitrary order of applying the
123	>	noncommuting rotation operators.\cite{Beard2003} Based on the
124	>	observation the momentum relaxation time is much less than the time
125	>	step, one may ignore the inertia in Brownian dynamics.  However, the
126	>	assumption of zero average acceleration is not always true for
127	>	cooperative motion which is common in proteins. An inertial Brownian
128	>	dynamics (IBD) was proposed to address this issue by adding an
129	>	inertial correction term.\cite{Beard2000} As a complement to IBD,
130	>	which has a lower bound in time step because of the inertial
131	>	relaxation time, long-time-step inertial dynamics (LTID) can be used
132	>	to investigate the inertial behavior of linked polymer segments in a
133	>	low friction regime.\cite{Beard2000} LTID can also deal with the
134	>	rotational dynamics for nonskew bodies without translation-rotation
135	>	coupling by separating the translation and rotation motion and taking
136	>	advantage of the analytical solution of hydrodynamic
137	>	properties. However, typical nonskew bodies like cylinders and
138	>	ellipsoids are inadequate to represent most complex macromolecular
139	>	assemblies. Therefore, the goal of this work is to adapt some of the
140	>	hydrodynamic methodologies developed to treat Brownian motion of
141	>	complex assemblies into a Langevin integrator for rigid bodies with
142	>	arbitrary shapes.
143
144	+	\subsection{Rigid Body Dynamics}
145	+	Rigid bodies are frequently involved in the modeling of large
146	+	collections of particles that move as a single unit.  In molecular
147	+	simulations, rigid bodies have been used to simplify protein-protein
148	+	docking,\cite{Gray2003} and lipid bilayer
149	+	simulations.\cite{SunX._jp0762020} Many of the water models in common
150	+	use are also rigid-body
151	+	models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are
152	+	typically evolved in molecular dynamics simulations using constraints
153	+	rather than rigid body equations of motion.
154	+
155	+	Euler angles are a natural choice to describe the rotational degrees
156	+	of freedom.  However, due to $\frac{1}{\sin \theta}$ singularities, the
157	+	numerical integration of corresponding equations of these motion can
158	+	become inaccurate (and inefficient).  Although the use of multiple
159	+	sets of Euler angles can overcome this problem,\cite{Barojas1973} the
160	+	computational penalty and the loss of angular momentum conservation
161	+	remain.  A singularity-free representation utilizing quaternions was
162	+	developed by Evans in 1977.\cite{Evans1977} The Evans quaternion
163	+	approach uses a nonseparable Hamiltonian, and this has prevented
164	+	symplectic algorithms from being utilized until very
165	+	recently.\cite{Miller2002}
166	+
167	+	Another approach is the application of holonomic constraints to the
168	+	atoms belonging to the rigid body.  Each atom moves independently
169	+	under the normal forces deriving from potential energy and constraints
170	+	are used to guarantee rigidity. However, due to their iterative
171	+	nature, the SHAKE and RATTLE algorithms converge very slowly when the
172	+	number of constraints (and the number of particles that belong to the
173	+	rigid body) increases.\cite{Ryckaert1977,Andersen1983}
174	+
175	+	In order to develop a stable and efficient integration scheme that
176	+	preserves most constants of the motion in microcanonical simulations,
177	+	symplectic propagators are necessary.  By introducing a conjugate
178	+	momentum to the rotation matrix ${\bf Q}$ and re-formulating
179	+	Hamilton's equations, a symplectic orientational integrator,
180	+	RSHAKE,\cite{Kol1997} was proposed to evolve rigid bodies on a
181	+	constraint manifold by iteratively satisfying the orthogonality
182	+	constraint ${\bf Q}^T {\bf Q} = 1$.  An alternative method using the
183	+	quaternion representation was developed by Omelyan.\cite{Omelyan1998}
184	+	However, both of these methods are iterative and suffer from some
185	+	related inefficiencies. A symplectic Lie-Poisson integrator for rigid
186	+	bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes
187	+	most of the limitations mentioned above and is therefore the basis for
188	+	our Langevin integrator.
189	+
190		The goal of the present work is to develop a Langevin dynamics
191	<	algorithm for arbitrary-shaped rigid particles by integrating the
192	<	accurate estimation of friction tensor from hydrodynamics theory
193	<	into the sophisticated rigid body dynamics algorithms.
191	>	algorithm for arbitrary-shaped rigid particles by integrating an
192	>	accurate estimate of the friction tensor from hydrodynamics theory
193	>	into a stable and efficient rigid body dynamics propagator.  In the
194	>	sections below, we review some of the theory of hydrodynamic tensors
195	>	developed primarily for Brownian simulations of multi-particle
196	>	systems, we then present our integration method for a set of
197	>	generalized Langevin equations of motion, and we compare the behavior
198	>	of the new Langevin integrator to dynamical quantities obtained via
199	>	explicit solvent molecular dynamics.
200
201	<	\subsection{\label{introSection:frictionTensor}Friction Tensor}
202	<	Theoretically, the friction kernel can be determined using the
203	<	velocity autocorrelation function. However, this approach becomes
204	<	impractical when the system becomes more and more complicated.
205	<	Instead, various approaches based on hydrodynamics have been
206	<	developed to calculate the friction coefficients. In general, the
207	<	friction tensor $\Xi$ is a $6\times 6$ matrix given by
208	<	\[
209	<	\Xi  = \left( {\begin{array}{*{20}c}
210	<	{\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
211	<	{\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
212	<	\end{array}} \right).
213	<	\]
214	<	Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$
215	<	translational friction tensor and rotational resistance (friction)
216	<	tensor respectively, while ${\Xi^{tr} }$ is translation-rotation
217	<	coupling tensor and $ {\Xi^{rt} }$ is rotation-translation coupling
218	<	tensor. When a particle moves in a fluid, it may experience friction
219	<	force or torque along the opposite direction of the velocity or
220	<	angular velocity,
221	<	\[
201	>	\subsection{\label{introSection:frictionTensor}The Friction Tensor}
202	>	Theoretically, a complete friction kernel for a solute particle can be
203	>	determined using the velocity autocorrelation function from a
204	>	simulation with explicit solvent molecules. However, this approach
205	>	becomes impractical when the solute becomes complex.  Instead, various
206	>	approaches based on hydrodynamics have been developed to calculate
207	>	static friction coefficients. In general, the friction tensor $\Xi$ is
208	>	a $6\times 6$ matrix given by
209	>	\begin{equation}
210	>	\Xi  = \left( \begin{array}{*{20}c}
211	>	\Xi^{tt} & \Xi^{rt}  \\
212	>	\Xi^{tr} & \Xi^{rr}  \\
213	>	\end{array} \right).
214	>	\end{equation}
215	>	Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and
216	>	rotational resistance (friction) tensors respectively, while
217	>	$\Xi^{tr}$ is translation-rotation coupling tensor and $\Xi^{rt}$ is
218	>	rotation-translation coupling tensor. When a particle moves in a
219	>	fluid, it may experience a friction force ($\mathbf{f}_f$) and torque
220	>	($\mathbf{\tau}_f$) in opposition to the velocity ($\mathbf{v}$) and
221	>	body-fixed angular velocity ($\mathbf{\omega}$),
222	>	\begin{equation}
223		\left( \begin{array}{l}
224	<	F_R  \\
225	<	\tau _R  \\
226	<	\end{array} \right) =  - \left( {\begin{array}{*{20}c}
227	<	{\Xi ^{tt} } & {\Xi ^{rt} }  \\
228	<	{\Xi ^{tr} } & {\Xi ^{rr} }  \\
229	<	\end{array}} \right)\left( \begin{array}{l}
230	<	v \\
231	<	w \\
232	<	\end{array} \right)
233	<	\]
234	<	where $F_r$ is the friction force and $\tau _R$ is the friction
235	<	torque.
224	>	\mathbf{f}_f  \\
225	>	\mathbf{\tau}_f  \\
226	>	\end{array} \right) =  - \left( \begin{array}{*{20}c}
227	>	\Xi^{tt} & \Xi^{rt}  \\
228	>	\Xi^{tr} & \Xi^{rr}  \\
229	>	\end{array} \right)\left( \begin{array}{l}
230	>	\mathbf{v} \\
231	>	\mathbf{\omega} \\
232	>	\end{array} \right).
233	>	\end{equation}
234	>	For an arbitrary body moving in a fluid, Peters has derived a set of
235	>	fluctuation-dissipation relations for the friction
236	>	tensors,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk}
237	>	\begin{eqnarray}
238	>	\Xi^{tt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
239	>	F}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}^2
240	>	\right] ds \\
241	>	\notag \\
242	>	\Xi^{tr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
243	>	F}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}
244	>	\langle {\bf \tau} \rangle_{eq} \right] ds \\
245	>	\notag \\
246	>	\Xi^{rt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
247	>	\tau}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}
248	>	\langle {\bf F} \rangle_{eq} \right] ds \\
249	>	\notag \\
250	>	\Xi^{rr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
251	>	\tau}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}^2
252	>	\right] ds
253	>	\end{eqnarray}
254	>	In these expressions, the forces (${\bf F}$) and torques (${\bf
255	>	\tau}$) are those that arise solely from the interactions of the body with
256	>	the surrounding fluid. For a single solute body in an isotropic fluid,
257	>	the average forces and torques in these expressions ($\langle {\bf F}
258	>	\rangle_{eq}$ and $\langle {\bf \tau} \rangle_{eq}$)
259	>	vanish, and one obtains the simpler force-torque correlation formulae
260	>	of Nienhuis.\cite{Nienhuis:1970lr} Molecular dynamics simulations with
261	>	explicit solvent molecules can be used to obtain estimates of the
262	>	friction tensors with these formulae. In practice, however, one needs
263	>	relatively long simulations with frequently-stored force and torque
264	>	information to compute friction tensors, and this becomes
265	>	prohibitively expensive when there are large numbers of large solute
266	>	particles.  For bodies with simple shapes, there are a number of
267	>	approximate expressions that allow computation of these tensors
268	>	without the need for expensive simulations that utilize explicit
269	>	solvent particles.
270
271		\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
272	<
273	<	For a spherical particle with slip boundary conditions, the
274	<	translational and rotational friction constant can be calculated
275	<	from Stoke's law,
276	<	\[
277	<	\Xi ^{tt}  = \left( {\begin{array}{*{20}c}
278	<	{6\pi \eta R} & 0 & 0  \\
279	<	0 & {6\pi \eta R} & 0  \\
280	<	0 & 0 & {6\pi \eta R}  \\
281	<	\end{array}} \right)
282	<	\]
272	>	For a spherical body under ``stick'' boundary conditions, the
273	>	translational and rotational friction tensors can be estimated from
274	>	Stokes' law,
275	>	\begin{equation}
276	>	\label{eq:StokesTranslation}
277	>	\Xi^{tt}  = \left( \begin{array}{*{20}c}
278	>	{6\pi \eta \rho} & 0 & 0  \\
279	>	0 & {6\pi \eta \rho} & 0  \\
280	>	0 & 0 & {6\pi \eta \rho}  \\
281	>	\end{array} \right)
282	>	\end{equation}
283		and
284	<	\[
285	<	\Xi ^{rr}  = \left( {\begin{array}{*{20}c}
286	<	{8\pi \eta R^3 } & 0 & 0  \\
287	<	0 & {8\pi \eta R^3 } & 0  \\
288	<	0 & 0 & {8\pi \eta R^3 }  \\
289	<	\end{array}} \right)
290	<	\]
291	<	where $\eta$ is the viscosity of the solvent and $R$ is the
292	<	hydrodynamic radius.
284	>	\begin{equation}
285	>	\label{eq:StokesRotation}
286	>	\Xi^{rr}  = \left( \begin{array}{*{20}c}
287	>	{8\pi \eta \rho^3 } & 0 & 0  \\
288	>	0 & {8\pi \eta \rho^3 } & 0  \\
289	>	0 & 0 & {8\pi \eta \rho^3 }  \\
290	>	\end{array} \right)
291	>	\end{equation}
292	>	where $\eta$ is the viscosity of the solvent and $\rho$ is the
293	>	hydrodynamic radius.  The presence of the rotational resistance tensor
294	>	implies that the spherical body has internal structure and
295	>	orientational degrees of freedom that must be propagated in time.  For
296	>	non-structured spherical bodies (i.e. the atoms in a traditional
297	>	molecular dynamics simulation) these degrees of freedom do not exist.
298
299		Other non-spherical shapes, such as cylinders and ellipsoids, are
300	<	widely used as references for developing new hydrodynamics theory,
300	>	widely used as references for developing new hydrodynamic theories,
301		because their properties can be calculated exactly. In 1936, Perrin
302	<	extended Stokes's law to general ellipsoids, also called a triaxial
303	<	ellipsoid, which is given in Cartesian coordinates
304	<	by\cite{Perrin1934, Perrin1936}
305	<	\[
306	<	\frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
307	<	}} = 1
308	<	\]
309	<	where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
310	<	due to the complexity of the elliptic integral, only the ellipsoid
311	<	with the restriction of two axes being equal, \textit{i.e.}
312	<	prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
313	<	exactly. Introducing an elliptic integral parameter $S$ for prolate
314	<	ellipsoids :
315	<	\[
316	<	S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
317	<	} }}{b},
318	<	\]
319	<	and oblate ellipsoids:
320	<	\[
321	<	S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
322	<	}}{a},
323	<	\]
324	<	one can write down the translational and rotational resistance
325	<	tensors
326	<	\begin{eqnarray*}
327	<	\Xi _a^{tt}  & = & 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}}. \\
328	<	\Xi _b^{tt}  & = & \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S +
329	<	2a}},
330	<	\end{eqnarray*}
331	<	and
332	<	\begin{eqnarray*}
266	<	\Xi _a^{rr} & = & \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}}, \\
267	<	\Xi _b^{rr} & = & \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}}.
268	<	\end{eqnarray*}
302	>	extended Stokes' law to general
303	>	ellipsoids,\cite{Perrin1934,Perrin1936} described in Cartesian
304	>	coordinates as
305	>	\begin{equation}
306	>	\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1.
307	>	\end{equation}
308	>	Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the
309	>	complexity of the elliptic integral, only uniaxial ellipsoids, either
310	>	prolate ($a \ge b = c$) or oblate ($a < b = c$), were solved
311	>	exactly. Introducing an elliptic integral parameter $S$ for prolate,
312	>	\begin{equation}
313	>	S = \frac{2}{\sqrt{a^2  - b^2}} \ln \frac{a + \sqrt{a^2  - b^2}}{b},
314	>	\end{equation}
315	>	and oblate,
316	>	\begin{equation}
317	>	S = \frac{2}{\sqrt {b^2  - a^2 }} \arctan \frac{\sqrt {b^2  - a^2}}{a},
318	>	\end{equation}
319	>	ellipsoids, it is possible to write down exact solutions for the
320	>	resistance tensors.  As is the case for spherical bodies, the translational,
321	>	\begin{eqnarray}
322	>	\Xi_a^{tt}  & = & 16\pi \eta \frac{a^2  - b^2}{(2a^2  - b^2 )S - 2a}. \\
323	>	\Xi_b^{tt} =  \Xi_c^{tt} & = & 32\pi \eta \frac{a^2  - b^2 }{(2a^2 - 3b^2 )S + 2a},
324	>	\end{eqnarray}
325	>	and rotational,
326	>	\begin{eqnarray}
327	>	\Xi_a^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^2  - b^2 )b^2}{2a - b^2 S}, \\
328	>	\Xi_b^{rr} = \Xi_c^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^4  - b^4)}{(2a^2  - b^2 )S - 2a}
329	>	\end{eqnarray}
330	>	resistance tensors are diagonal $3 \times 3$ matrices. For both
331	>	spherical and ellipsoidal particles, the translation-rotation and
332	>	rotation-translation coupling tensors are zero.
333
334		\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
335	+	Other than the fluctuation dissipation formulae given by
336	+	Peters,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk} there are no
337	+	analytic solutions for the friction tensor for rigid molecules of
338	+	arbitrary shape. The ellipsoid of revolution and general triaxial
339	+	ellipsoid models have been widely used to approximate the hydrodynamic
340	+	properties of rigid bodies. However, the mapping from all possible
341	+	ellipsoidal spaces ($r$-space) to all possible combinations of
342	+	rotational diffusion coefficients ($D$-space) is not
343	+	unique.\cite{Wegener1979} Additionally, because there is intrinsic
344	+	coupling between translational and rotational motion of {\it skew}
345	+	rigid bodies, general ellipsoids are not always suitable for modeling
346	+	rigid molecules.  A number of studies have been devoted to determining
347	+	the friction tensor for irregular shapes using methods in which the
348	+	molecule of interest is modeled with a combination of
349	+	spheres\cite{Carrasco1999} and the hydrodynamic properties of the
350	+	molecule are then calculated using a set of two-point interaction
351	+	tensors.  We have found the {\it bead} and {\it rough shell} models of
352	+	Carrasco and Garc\'{i}a de la Torre to be the most useful of these
353	+	methods,\cite{Carrasco1999} and we review the basic outline of the
354	+	rough shell approach here.  A more thorough explanation can be found
355	+	in Ref. \citen{Carrasco1999}.
356
357	<	Unlike spherical and other simply shaped molecules, there is no
358	<	analytical solution for the friction tensor for arbitrarily shaped
359	<	rigid molecules. The ellipsoid of revolution model and general
360	<	triaxial ellipsoid model have been used to approximate the
361	<	hydrodynamic properties of rigid bodies. However, since the mapping
277	<	from all possible ellipsoidal spaces, $r$-space, to all possible
278	<	combination of rotational diffusion coefficients, $D$-space, is not
279	<	unique\cite{Wegener1979} as well as the intrinsic coupling between
280	<	translational and rotational motion of rigid bodies, general
281	<	ellipsoids are not always suitable for modeling arbitrarily shaped
282	<	rigid molecules. A number of studies have been devoted to
283	<	determining the friction tensor for irregularly shaped rigid bodies
284	<	using more advanced methods where the molecule of interest was
285	<	modeled by a combinations of spheres\cite{Carrasco1999} and the
286	<	hydrodynamics properties of the molecule can be calculated using the
287	<	hydrodynamic interaction tensor. Let us consider a rigid assembly of
288	<	$N$ beads immersed in a continuous medium. Due to hydrodynamic
289	<	interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
290	<	than its unperturbed velocity $v_i$,
291	<	\[
292	<	v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
293	<	\]
294	<	where $F_i$ is the frictional force, and $T_{ij}$ is the
295	<	hydrodynamic interaction tensor. The friction force of $i$th bead is
296	<	proportional to its ``net'' velocity
357	>	Consider a rigid assembly of $N$ small beads moving through a
358	>	continuous medium.  Due to hydrodynamic interactions between the
359	>	beads, the net velocity of the $i^\mathrm{th}$ bead relative to the
360	>	medium, ${\bf v}'_i$, is different than its unperturbed velocity ${\bf
361	>	v}_i$,
362		\begin{equation}
363	<	F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
299	<	\label{introEquation:tensorExpression}
363	>	{\bf v}'_i  = {\bf v}_i  - \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }
364		\end{equation}
365	<	This equation is the basis for deriving the hydrodynamic tensor. In
366	<	1930, Oseen and Burgers gave a simple solution to
367	<	Eq.~\ref{introEquation:tensorExpression}
365	>	where ${\bf F}_j$ is the frictional force on the medium due to bead $j$, and
366	>	${\bf T}_{ij}$ is the hydrodynamic interaction tensor between the two beads.
367	>	The frictional force felt by the $i^\mathrm{th}$ bead is proportional to
368	>	its net velocity
369		\begin{equation}
370	<	T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
371	<	R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
370	>	{\bf F}_i  = \xi_i {\bf v}_i  - \xi_i \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }.
371	>	\label{introEquation:tensorExpression}
372		\end{equation}
373	<	Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
374	<	A second order expression for element of different size was
375	<	introduced by Rotne and Prager\cite{Rotne1969} and improved by
373	>	Eq. (\ref{introEquation:tensorExpression}) defines the two-point
374	>	hydrodynamic tensor, ${\bf T}_{ij}$.  There have been many proposed
375	>	solutions to this equation, including the simple solution given by
376	>	Oseen and Burgers in 1930 for two beads of identical radius.  A second
377	>	order expression for beads of different hydrodynamic radii was
378	>	introduced by Rotne and Prager,\cite{Rotne1969} and improved by
379		Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977}
380		\begin{equation}
381	<	T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
382	<	\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
383	<	_i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
384	<	\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
381	>	{\bf T}_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {{\bf I} +
382	>	\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\rho
383	>	_i^2  + \rho_j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} -
384	>	\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
385		\label{introEquation:RPTensorNonOverlapped}
386		\end{equation}
387	<	Both of the Eq.~\ref{introEquation:oseenTensor} and
388	<	Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption
389	<	$R_{ij} \ge \sigma _i  + \sigma _j$. An alternative expression for
390	<	overlapping beads with the same radius, $\sigma$, is given by
387	>	Here ${\bf R}_{ij}$ is the distance vector between beads $i$ and $j$.  Both
388	>	the Oseen-Burgers tensor and
389	>	Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption that
390	>	the beads do not overlap ($R_{ij} \ge \rho_i + \rho_j$).
391	>
392	>	To calculate the resistance tensor for a body represented as the union
393	>	of many non-overlapping beads, we first pick an arbitrary origin $O$
394	>	and then construct a $3N \times 3N$ supermatrix consisting of $N
395	>	\times N$ ${\bf B}_{ij}$ blocks
396		\begin{equation}
397	<	T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
398	<	\frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
399	<	\frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
400	<	\label{introEquation:RPTensorOverlapped}
397	>	{\bf B} = \left( \begin{array}{*{20}c}
398	>	{\bf B}_{11} &  \ldots  & {\bf B}_{1N}   \\
399	>	\vdots  &  \ddots  &  \vdots   \\
400	>	{\bf B}_{N1} &  \cdots  & {\bf B}_{NN}
401	>	\end{array} \right)
402		\end{equation}
403	<	To calculate the resistance tensor at an arbitrary origin $O$, we
404	<	construct a $3N \times 3N$ matrix consisting of $N \times N$
331	<	$B_{ij}$ blocks
403	>	${\bf B}_{ij}$ is a version of the hydrodynamic tensor which includes the
404	>	self-contributions for spheres,
405		\begin{equation}
406	<	B = \left( {\begin{array}{*{20}c}
407	<	{B_{11} } &  \ldots  & {B_{1N} }  \\
335	<	\vdots  &  \ddots  &  \vdots   \\
336	<	{B_{N1} } &  \cdots  & {B_{NN} }  \\
337	<	\end{array}} \right),
406	>	{\bf B}_{ij}  = \delta _{ij} \frac{{\bf I}}{{6\pi \eta R_{ij}}} + (1 - \delta_{ij}
407	>	){\bf T}_{ij}
408		\end{equation}
409	<	where $B_{ij}$ is given by
410	<	\[
411	<	B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
412	<	)T_{ij}
413	<	\]
414	<	where $\delta _{ij}$ is the Kronecker delta function. Inverting the
415	<	$B$ matrix, we obtain
416	<	\[
417	<	C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
418	<	{C_{11} } &  \ldots  & {C_{1N} }  \\
419	<	\vdots  &  \ddots  &  \vdots   \\
420	<	{C_{N1} } &  \cdots  & {C_{NN} }  \\
421	<	\end{array}} \right),
422	<	\]
423	<	which can be partitioned into $N \times N$ $3 \times 3$ block
424	<	$C_{ij}$. With the help of $C_{ij}$ and the skew matrix $U_i$
425	<	\[
426	<	U_i  = \left( {\begin{array}{*{20}c}
427	<	0 & { - z_i } & {y_i }  \\
428	<	{z_i } & 0 & { - x_i }  \\
359	<	{ - y_i } & {x_i } & 0  \\
360	<	\end{array}} \right)
361	<	\]
409	>	where $\delta_{ij}$ is the Kronecker delta function. Inverting the
410	>	${\bf B}$ matrix, we obtain
411	>	\begin{equation}
412	>	{\bf C} = {\bf B}^{ - 1}  = \left(\begin{array}{*{20}c}
413	>	{\bf C}_{11} &  \ldots  & {\bf C}_{1N} \\
414	>	\vdots  &  \ddots  &  \vdots   \\
415	>	{\bf C}_{N1} &  \cdots  & {\bf C}_{NN}
416	>	\end{array} \right),
417	>	\end{equation}
418	>	which can be partitioned into $N \times N$ blocks labeled ${\bf C}_{ij}$.
419	>	(Each of the ${\bf C}_{ij}$ blocks is a $3 \times 3$ matrix.)  Using the
420	>	skew matrix,
421	>	\begin{equation}
422	>	{\bf U}_i  = \left(\begin{array}{*{20}c}
423	>	0 & -z_i & y_i \\
424	>	z_i &  0   & - x_i \\
425	>	-y_i & x_i & 0
426	>	\end{array}\right)
427	>	\label{eq:skewMatrix}
428	>	\end{equation}
429		where $x_i$, $y_i$, $z_i$ are the components of the vector joining
430	<	bead $i$ and origin $O$, the elements of resistance tensor at
431	<	arbitrary origin $O$ can be written as
430	>	bead $i$ and origin $O$, the elements of the resistance tensor (at the
431	>	arbitrary origin $O$) can be written as
432		\begin{eqnarray}
366	–	\Xi _{}^{tt}  & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
367	–	\Xi _{}^{tr}  & = & \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
368	–	\Xi _{}^{rr}  & = &  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } }
369	–	U_j  + 6 \eta V {\bf I}. \notag
433		\label{introEquation:ResistanceTensorArbitraryOrigin}
434	+	\Xi^{tt}  & = & \sum\limits_i {\sum\limits_j {{\bf C}_{ij} } } \notag , \\
435	+	\Xi^{tr}  = \Xi _{}^{rt}  & = & \sum\limits_i {\sum\limits_j {{\bf U}_i {\bf C}_{ij} } } , \\
436	+	\Xi^{rr}  & = &  -\sum\limits_i \sum\limits_j {\bf U}_i {\bf C}_{ij} {\bf U}_j + 6 \eta V {\bf I}. \notag
437		\end{eqnarray}
438	<	The final term in the expression for $\Xi^{rr}$ is correction that
439	<	accounts for errors in the rotational motion of certain kinds of bead
440	<	models. The additive correction uses the solvent viscosity ($\eta$)
441	<	as well as the total volume of the beads that contribute to the
376	<	hydrodynamic model,
438	>	The final term in the expression for $\Xi^{rr}$ is a correction that
439	>	accounts for errors in the rotational motion of the bead models. The
440	>	additive correction uses the solvent viscosity ($\eta$) as well as the
441	>	total volume of the beads that contribute to the hydrodynamic model,
442		\begin{equation}
443	<	V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3,
443	>	V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3,
444		\end{equation}
445	<	where $\sigma_i$ is the radius of bead $i$.  This correction term was
445	>	where $\rho_i$ is the radius of bead $i$.  This correction term was
446		rigorously tested and compared with the analytical results for
447	<	two-sphere and ellipsoidal systems by Garcia de la Torre and
447	>	two-sphere and ellipsoidal systems by Garc\'{i}a de la Torre and
448		Rodes.\cite{Torre:1983lr}
449
450	<
451	<	The resistance tensor depends on the origin to which they refer. The
452	<	proper location for applying the friction force is the center of
453	<	resistance (or center of reaction), at which the trace of rotational
454	<	resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
455	<	Mathematically, the center of resistance is defined as an unique
456	<	point of the rigid body at which the translation-rotation coupling
392	<	tensors are symmetric,
450	>	In general, resistance tensors depend on the origin at which they were
451	>	computed.  However, the proper location for applying the friction
452	>	force is the center of resistance, the special point at which the
453	>	trace of rotational resistance tensor, $\Xi^{rr}$ reaches a minimum
454	>	value.  Mathematically, the center of resistance can also be defined
455	>	as the unique point for a rigid body at which the translation-rotation
456	>	coupling tensors are symmetric,
457		\begin{equation}
458	<	\Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
458	>	\Xi^{tr}  = \left(\Xi^{tr}\right)^T
459		\label{introEquation:definitionCR}
460		\end{equation}
461	<	From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
462	<	we can easily derive that the translational resistance tensor is
463	<	origin independent, while the rotational resistance tensor and
461	>	From Eq. \ref{introEquation:ResistanceTensorArbitraryOrigin}, we can
462	>	easily derive that the {\it translational} resistance tensor is origin
463	>	independent, while the rotational resistance tensor and
464		translation-rotation coupling resistance tensor depend on the
465	<	origin. Given the resistance tensor at an arbitrary origin $O$, and
466	<	a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
467	<	obtain the resistance tensor at $P$ by
468	<	\begin{equation}
469	<	\begin{array}{l}
470	<	\Xi _P^{tt}  = \Xi _O^{tt}  \\
471	<	\Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
472	<	\Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{{tr} ^{^T }}  \\
473	<	\end{array}
474	<	\label{introEquation:resistanceTensorTransformation}
475	<	\end{equation}
476	<	where
477	<	\[
478	<	U_{OP}  = \left( {\begin{array}{*{20}c}
479	<	0 & { - z_{OP} } & {y_{OP} }  \\
480	<	{z_i } & 0 & { - x_{OP} }  \\
481	<	{ - y_{OP} } & {x_{OP} } & 0  \\
482	<	\end{array}} \right)
483	<	\]
484	<	Using Eq.~\ref{introEquation:definitionCR} and
485	<	Eq.~\ref{introEquation:resistanceTensorTransformation}, one can
486	<	locate the position of center of resistance,
487	<	\begin{eqnarray*}
488	<	\left( \begin{array}{l}
489	<	x_{OR}  \\
490	<	y_{OR}  \\
491	<	z_{OR}  \\
492	<	\end{array} \right) & = &\left( {\begin{array}{*{20}c}
493	<	{(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
494	<	{ - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
495	<	{ - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
496	<	\end{array}} \right)^{ - 1}  \\
497	<	& & \left( \begin{array}{l}
434	<	(\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
435	<	(\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
436	<	(\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
437	<	\end{array} \right) \\
438	<	\end{eqnarray*}
439	<	where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
465	>	origin. Given the resistance tensor at an arbitrary origin $O$, and a
466	>	vector ,${\bf r}_{OP} = (x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we
467	>	can obtain the resistance tensor at $P$ by
468	>	\begin{eqnarray}
469	>	\label{introEquation:resistanceTensorTransformation}
470	>	\Xi_P^{tt}  & = & \Xi_O^{tt}  \notag \\
471	>	\Xi_P^{tr}  = \Xi_P^{rt}  & = & \Xi_O^{tr}  - {\bf U}_{OP} \Xi _O^{tt}  \\
472	>	\Xi_P^{rr}  & = &\Xi_O^{rr}  - {\bf U}_{OP} \Xi_O^{tt} {\bf U}_{OP}
473	>	+ \Xi_O^{tr} {\bf U}_{OP}  - {\bf U}_{OP} \left( \Xi_O^{tr}
474	>	\right)^{^T} \notag
475	>	\end{eqnarray}
476	>	where ${\bf U}_{OP}$ is the skew matrix (Eq. (\ref{eq:skewMatrix}))
477	>	for the vector between the origin $O$ and the point $P$. Using
478	>	Eqs.~\ref{introEquation:definitionCR}~and~\ref{introEquation:resistanceTensorTransformation},
479	>	one can locate the position of center of resistance,
480	>	\begin{equation*}
481	>	\left(\begin{array}{l}
482	>	x_{OR} \\
483	>	y_{OR} \\
484	>	z_{OR}
485	>	\end{array}\right) =
486	>	\left(\begin{array}{*{20}c}
487	>	(\Xi_O^{rr})_{yy} + (\Xi_O^{rr})_{zz} & -(\Xi_O^{rr})_{xy} & -(\Xi_O^{rr})_{xz} \\
488	>	-(\Xi_O^{rr})_{xy} & (\Xi_O^{rr})_{zz} + (\Xi_O^{rr})_{xx} & -(\Xi_O^{rr})_{yz} \\
489	>	-(\Xi_O^{rr})_{xz} & -(\Xi_O^{rr})_{yz} & (\Xi_O^{rr})_{xx} + (\Xi_O^{rr})_{yy} \\
490	>	\end{array}\right)^{-1}
491	>	\left(\begin{array}{l}
492	>	(\Xi_O^{tr})_{yz} - (\Xi_O^{tr})_{zy} \\
493	>	(\Xi_O^{tr})_{zx} - (\Xi_O^{tr})_{xz} \\
494	>	(\Xi_O^{tr})_{xy} - (\Xi_O^{tr})_{yx}
495	>	\end{array}\right)
496	>	\end{equation*}
497	>	where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector
498		joining center of resistance $R$ and origin $O$.
499
500	+	For a general rigid molecular substructure, finding the $6 \times 6$
501	+	resistance tensor can be a computationally demanding task.  First, a
502	+	lattice of small beads that extends well beyond the boundaries of the
503	+	rigid substructure is created.  The lattice is typically composed of
504	+	0.25 \AA\ beads on a dense FCC lattice.  The lattice constant is taken
505	+	to be the bead diameter, so that adjacent beads are touching, but do
506	+	not overlap. To make a shape corresponding to the rigid structure,
507	+	beads that sit on lattice sites that are outside the van der Waals
508	+	radii of all of the atoms comprising the rigid body are excluded from
509	+	the calculation.
510
511	+	For large structures, most of the beads will be deep within the rigid
512	+	body and will not contribute to the hydrodynamic tensor.  In the {\it
513	+	rough shell} approach, beads which have all of their lattice neighbors
514	+	inside the structure are considered interior beads, and are removed
515	+	from the calculation.  After following this procedure, only those
516	+	beads in direct contact with the van der Waals surface of the rigid
517	+	body are retained.  For reasonably large molecular structures, this
518	+	truncation can still produce bead assemblies with thousands of
519	+	members.
520	+
521	+	If all of the {\it atoms} comprising the rigid substructure are
522	+	spherical and non-overlapping, the tensor in
523	+	Eq.~(\ref{introEquation:RPTensorNonOverlapped}) may be used directly
524	+	using the atoms themselves as the hydrodynamic beads.  This is a
525	+	variant of the {\it bead model} approach of Carrasco and Garc\'{i}a de
526	+	la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$
527	+	matrix can be quite small, and the calculation of the hydrodynamic
528	+	tensor is straightforward.
529	+
530	+	In general, the inversion of the ${\bf B}$ matrix is the most
531	+	computationally demanding task.  This inversion is done only once for
532	+	each type of rigid structure.  We have used straightforward
533	+	LU-decomposition to solve the linear system and to obtain the elements
534	+	of ${\bf C}$. Once ${\bf C}$ has been obtained, the location of the
535	+	center of resistance ($R$) is found and the resistance tensor at this
536	+	point is calculated.  The $3 \times 1$ vector giving the location of
537	+	the rigid body's center of resistance and the $6 \times 6$ resistance
538	+	tensor are then stored for use in the Langevin dynamics calculation.
539	+	These quantities depend on solvent viscosity and temperature and must
540	+	be recomputed if different simulation conditions are required.
541	+
542		\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
543	+
544		Consider the Langevin equations of motion in generalized coordinates
545		\begin{equation}
546	<	M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
546	>	{\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) +
547	>	{\bf F}_{f}(t)  + {\bf F}_{r}(t)
548		\label{LDGeneralizedForm}
549		\end{equation}
550	<	where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
551	<	and moment of inertial) matrix and $V_i$ is a generalized velocity,
552	<	$V_i = V_i(v_i,\omega _i)$. The right side of
553	<	Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in
554	<	lab-fixed frame, systematic force $F_{s,i}$, dissipative force
555	<	$F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
556	<	system in Newtownian mechanics typically refers to lab-fixed frame,
557	<	it is also convenient to handle the rotation of rigid body in
558	<	body-fixed frame. Thus the friction and random forces are calculated
559	<	in body-fixed frame and converted back to lab-fixed frame by:
560	<	\[
561	<	\begin{array}{l}
562	<	F_{f,i}^l (t) = Q^T F_{f,i}^b (t), \\
563	<	F_{r,i}^l (t) = Q^T F_{r,i}^b (t). \\
463	<	\end{array}
464	<	\]
465	<	Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
466	<	the body-fixed velocity at center of resistance $v_{R,i}^b$ and
467	<	angular velocity $\omega _i$
550	>	where ${\bf M}$ is a $6 \times 6$ diagonal mass matrix (which
551	>	includes the mass of the rigid body as well as the moments of inertia
552	>	in the body-fixed frame) and ${\bf V}$ is a generalized velocity,
553	>	${\bf V} =
554	>	\left\{{\bf v},{\bf \omega}\right\}$. The right side of
555	>	Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a
556	>	system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf
557	>	F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution
558	>	of the system in Newtonian mechanics is typically done in the lab
559	>	frame, it is convenient to handle the dynamics of rigid bodies in
560	>	body-fixed frames. Thus the friction and random forces on each
561	>	substructure are calculated in a body-fixed frame and may converted
562	>	back to the lab frame using that substructure's rotation matrix (${\bf
563	>	Q}$):
564		\begin{equation}
565	<	F_{r,i}^b (t) = \left( \begin{array}{l}
566	<	f_{r,i}^b (t) \\
567	<	\tau _{r,i}^b (t) \\
568	<	\end{array} \right) =  - \left( {\begin{array}{*{20}c}
569	<	{\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
570	<	{\Xi _{R,c} } & {\Xi _{R,r} }  \\
571	<	\end{array}} \right)\left( \begin{array}{l}
572	<	v_{R,i}^b (t) \\
573	<	\omega _i (t) \\
565	>	{\bf F}_{f,r} =
566	>	\left( \begin{array}{c}
567	>	{\bf f}_{f,r} \\
568	>	{\bf \tau}_{f,r}
569	>	\end{array} \right)
570	>	=
571	>	\left( \begin{array}{c}
572	>	{\bf Q}^{T} {\bf f}^{~b}_{f,r} \\
573	>	{\bf Q}^{T} {\bf \tau}^{~b}_{f,r}
574	>	\end{array} \right)
575	>	\end{equation}
576	>	The body-fixed friction force, ${\bf F}_{f}^{~b}$, is proportional to
577	>	the (body-fixed) velocity at the center of resistance
578	>	${\bf v}_{R}^{~b}$ and the angular velocity ${\bf \omega}$
579	>	\begin{equation}
580	>	{\bf F}_{f}^{~b}(t) = \left( \begin{array}{l}
581	>	{\bf f}_{f}^{~b}(t) \\
582	>	{\bf \tau}_{f}^{~b}(t) \\
583	>	\end{array} \right) =  - \left( \begin{array}{*{20}c}
584	>	\Xi_{R}^{tt} & \Xi_{R}^{rt} \\
585	>	\Xi_{R}^{tr} & \Xi_{R}^{rr}    \\
586	>	\end{array} \right)\left( \begin{array}{l}
587	>	{\bf v}_{R}^{~b}(t) \\
588	>	{\bf \omega}(t) \\
589		\end{array} \right),
590		\end{equation}
591	<	while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
592	<	with zero mean and variance
591	>	while the random force, ${\bf F}_{r}$, is a Gaussian stochastic
592	>	variable with zero mean and variance,
593		\begin{equation}
594	<	\left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
595	<	\left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
596	<	2k_B T\Xi _R \delta (t - t'). \label{randomForce}
594	>	\left\langle {{\bf F}_{r}(t) ({\bf F}_{r}(t'))^T } \right\rangle  =
595	>	\left\langle {{\bf F}_{r}^{~b} (t) ({\bf F}_{r}^{~b} (t'))^T } \right\rangle  =
596	>	2 k_B T \Xi_R \delta(t - t'). \label{eq:randomForce}
597		\end{equation}
598	<	The equation of motion for $v_i$ can be written as
598	>	$\Xi_R$ is the $6\times6$ resistance tensor at the center of
599	>	resistance.  Once this tensor is known for a given rigid body (as
600	>	described in the previous section) obtaining a stochastic vector that
601	>	has the properties in Eq. (\ref{eq:randomForce}) can be done
602	>	efficiently by carrying out a one-time Cholesky decomposition to
603	>	obtain the square root matrix of the resistance tensor,
604	>	\begin{equation}
605	>	\Xi_R = {\bf S} {\bf S}^{T},
606	>	\label{eq:Cholesky}
607	>	\end{equation}
608	>	where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
609	>	vector with the statistics required for the random force can then be
610	>	obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which
611	>	has elements chosen from a Gaussian distribution, such that:
612		\begin{equation}
613	<	m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
614	<	f_{r,i}^l (t)
613	>	\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
614	>	{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
615		\end{equation}
616	<	Since the frictional force is applied at the center of resistance
617	<	which generally does not coincide with the center of mass, an extra
618	<	torque is exerted at the center of mass. Thus, the net body-fixed
619	<	frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
620	<	given by
616	>	where $\delta t$ is the timestep in use during the simulation. The
617	>	random force, ${\bf F}_{r}^{~b} = {\bf S} {\bf Z}$, can be shown to have the
618	>	correct properties required by Eq. (\ref{eq:randomForce}).
619	>
620	>	The equation of motion for the translational velocity of the center of
621	>	mass (${\bf v}$) can be written as
622		\begin{equation}
623	<	\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
623	>	m \dot{{\bf v}} (t) =  {\bf f}_{s}(t) + {\bf f}_{f}(t) +
624	>	{\bf f}_{r}(t)
625		\end{equation}
626	<	where $r_{MR}$ is the vector from the center of mass to the center
627	<	of the resistance. Instead of integrating the angular velocity in
628	<	lab-fixed frame, we consider the equation of angular momentum in
629	<	body-fixed frame
626	>	Since the frictional and random forces are applied at the center of
627	>	resistance, which generally does not coincide with the center of mass,
628	>	extra torques are exerted at the center of mass. Thus, the net
629	>	body-fixed torque at the center of mass, $\tau^{~b}(t)$,
630	>	is given by
631		\begin{equation}
632	<	\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t)
506	<	+ \tau _{r,i}^b(t)
632	>	\tau^{~b} \leftarrow \tau_{s}^{~b} + \tau_{f}^{~b} + \tau_{r}^{~b} + {\bf r}_{MR} \times \left( {\bf f}_{f}^{~b} + {\bf f}_{r}^{~b} \right)
633		\end{equation}
634	<	Embedding the friction terms into force and torque, one can
635	<	integrate the langevin equations of motion for rigid body of
636	<	arbitrary shape in a velocity-Verlet style 2-part algorithm, where
637	<	$h= \delta t$:
634	>	where ${\bf r}_{MR}$ is the vector from the center of mass to the center of
635	>	resistance. Instead of integrating the angular velocity in lab-fixed
636	>	frame, we consider the equation of motion for the angular momentum
637	>	(${\bf j}$) in the body-fixed frame
638	>	\begin{equation}
639	>	\frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t)
640	>	\end{equation}
641	>	Embedding the friction and random forces into the the total force and
642	>	torque, one can integrate the Langevin equations of motion for a rigid
643	>	body of arbitrary shape in a velocity-Verlet style 2-part algorithm,
644	>	where $h = \delta t$:
645
646	<	{\tt moveA:}
646	>	{\tt move A:}
647		\begin{align*}
648		{\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t)
649		+ \frac{h}{2} \left( {\bf f}(t) / m \right), \\
#	Line 519 | Line 652 | $h= \delta t$:
652		+ h  {\bf v}\left(t + h / 2 \right), \\
653		%
654		{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
655	<	+ \frac{h}{2} {\bf \tau}^b(t), \\
655	>	+ \frac{h}{2} {\bf \tau}^{~b}(t), \\
656		%
657	<	\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
657	>	{\bf Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
658		(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
659		\end{align*}
660	<	In this context, the $\mathrm{rotate}$ function is the reversible
661	<	product of the three body-fixed rotations,
660	>	In this context, $\overleftrightarrow{\mathsf{I}}$ is the diagonal
661	>	moment of inertia tensor, and the $\mathrm{rotate}$ function is the
662	>	reversible product of the three body-fixed rotations,
663		\begin{equation}
664		\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
665		\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
666		/ 2) \cdot \mathsf{G}_x(a_x /2),
667		\end{equation}
668		where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
669	<	rotates both the rotation matrix ($\mathsf{Q}$) and the body-fixed
669	>	rotates both the rotation matrix ($\mathbf{Q}$) and the body-fixed
670		angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
671		axis $\alpha$,
672		\begin{equation}
673		\mathsf{G}_\alpha( \theta ) = \left\{
674		\begin{array}{lcl}
675	<	\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
675	>	\mathbf{Q}(t) & \leftarrow & \mathbf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
676		{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
677		j}(0).
678		\end{array}
#	Line 560 | Line 694 | calculated at the new positions and orientations
694		\end{equation}
695		All other rotations follow in a straightforward manner. After the
696		first part of the propagation, the forces and body-fixed torques are
697	<	calculated at the new positions and orientations
697	>	calculated at the new positions and orientations.  The system forces
698	>	and torques are derivatives of the total potential energy function
699	>	($U$) with respect to the rigid body positions (${\bf r}$) and the
700	>	columns of the transposed rotation matrix ${\bf Q}^T = \left({\bf
701	>	u}_x, {\bf u}_y, {\bf u}_z \right)$:
702
703	<	{\tt doForces:}
703	>	{\tt Forces:}
704		\begin{align*}
705	<	{\bf f}(t + h) &\leftarrow
706	<	- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
705	>	{\bf f}_{s}(t + h) & \leftarrow
706	>	- \left(\frac{\partial U}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\
707		%
708	<	{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
709	<	\times \frac{\partial V}{\partial {\bf u}}, \\
708	>	{\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h)
709	>	\times \frac{\partial U}{\partial {\bf u}} \\
710		%
711	<	{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
712	<	\cdot {\bf \tau}^s(t + h).
711	>	{\bf v}^{b}_{R}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \left({\bf v}(t+h) + {\bf \omega}(t+h) \times {\bf r}_{MR} \right) \\
712	>	%
713	>	{\bf f}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tt} \cdot
714	>	{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rt} \cdot {\bf \omega}(t+h) \\
715	>	%
716	>	{\bf \tau}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tr} \cdot
717	>	{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rr} \cdot {\bf \omega}(t+h) \\
718	>	%
719	>	Z & \leftarrow  {\tt GaussianNormal}(2 k_B T / h, 6) \\
720	>	{\bf F}_{R,r}^{b}(t+h) & \leftarrow {\bf S} \cdot Z \\
721	>	%
722	>	{\bf f}(t+h) & \leftarrow {\bf f}_{s}(t+h) + \mathbf{Q}^{T}(t+h)
723	>	\cdot \left({\bf f}_{R,f}^{~b} + {\bf f}_{R,r}^{~b} \right) \\
724	>	%
725	>	\tau(t+h) & \leftarrow \tau_{s}(t+h) + \mathbf{Q}^{T}(t+h) \cdot \left(\tau_{R,f}^{~b} + \tau_{R,r}^{~b} \right) + {\bf r}_{MR} \times \left({\bf f}_{f}(t+h) + {\bf f}_{r}(t+h) \right) \\
726	>	\tau^{~b}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \tau(t+h) \\
727		\end{align*}
728	+	Frictional (and random) forces and torques must be computed at the
729	+	center of resistance, so there are additional steps required to find
730	+	the body-fixed velocity (${\bf v}_{R}^{~b}$) at this location.  Mapping
731	+	the frictional and random forces at the center of resistance back to
732	+	the center of mass also introduces an additional term in the torque
733	+	one obtains at the center of mass.
734	+
735		Once the forces and torques have been obtained at the new time step,
736		the velocities can be advanced to the same time value.
737
738	<	{\tt moveB:}
738	>	{\tt move B:}
739		\begin{align*}
740		{\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
741		\right)
#	Line 584 | Line 743 | the velocities can be advanced to the same time value.
743		%
744		{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2
745		\right)
746	<	+ \frac{h}{2} {\bf \tau}^b(t + h) .
746	>	+ \frac{h}{2} {\bf \tau}^{~b}(t + h) .
747		\end{align*}
748
749		\section{Validating the Method\label{sec:validating}}
#	Line 647 | Line 806 | simulation box.  All simulations were equilibrated usi
806		We performed reference microcanonical simulations with explicit
807		solvents for each of the different model system.  In each case there
808		was one solute model and 1929 solvent molecules present in the
809	<	simulation box.  All simulations were equilibrated using a
809	>	simulation box.  All simulations were equilibrated for 5 ns using a
810		constant-pressure and temperature integrator with target values of 300
811		K for the temperature and 1 atm for pressure.  Following this stage,
812	<	further equilibration and sampling was done in a microcanonical
813	<	ensemble.  Since the model bodies are typically quite massive, we were
814	<	able to use a time step of 25 fs.
812	>	further equilibration (5 ns) and sampling (10 ns) was done in a
813	>	microcanonical ensemble.  Since the model bodies are typically quite
814	>	massive, we were able to use a time step of 25 fs.
815
816		The model systems studied used both Lennard-Jones spheres as well as
817		uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne
818		potential was a single site model for the interactions of rigid
819	<	ellipsoidal molecules.\cite{Gay81} It can be thought of as a
819	>	ellipsoidal molecules.\cite{Gay1981} It can be thought of as a
820		modification of the Gaussian overlap model originally described by
821		Berne and Pechukas.\cite{Berne72} The potential is constructed in the
822		familiar form of the Lennard-Jones function using
823		orientation-dependent $\sigma$ and $\epsilon$ parameters,
824		\begin{equation*}
825	<	V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
826	<	r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
827	<	{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u
825	>	V_{ij}({{\bf \hat u}_i}, {{\bf \hat u}_j}, {{\bf \hat
826	>	r}_{ij}}) = 4\epsilon ({{\bf \hat u}_i}, {{\bf \hat u}_j},
827	>	{{\bf \hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u
828		}_i},
829	<	{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
830	<	-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
831	<	{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
829	>	{{\bf \hat u}_j}, {{\bf \hat r}_{ij}})+\sigma_0}\right)^{12}
830	>	-\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u}_i}, {{\bf \hat u}_j},
831	>	{{\bf \hat r}_{ij}})+\sigma_0}\right)^6\right]
832		\label{eq:gb}
833		\end{equation*}
834
#	Line 686 | Line 845 | are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGe
845		Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e /
846		\epsilon^s$, describes the ratio between the well depths in the {\it
847		end-to-end} and side-by-side configurations.  Details of the potential
848	<	are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an
848	>	are given elsewhere,\cite{Luckhurst90,Golubkov06,SunX._jp0762020} and an
849		excellent overview of the computational methods that can be used to
850		efficiently compute forces and torques for this potential can be found
851		in Ref. \citen{Golubkov06}
#	Line 721 | Line 880 | A similar form exists for the bulk viscosity
880		\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}.
881		\label{eq:shear}
882		\end{equation}
883	<	A similar form exists for the bulk viscosity
884	<	\begin{equation}
726	<	\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left(
727	<	\int_{t_0}^{t_0 + t}
728	<	\left(P\left(t'\right)-\left\langle P \right\rangle \right)dt'
729	<	\right)^2 \right\rangle_{t_0}.
730	<	\end{equation}
731	<	Alternatively, the shear viscosity can also be calculated using a
732	<	Green-Kubo formula with the off-diagonal pressure tensor correlation function,
733	<	\begin{equation}
734	<	\eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0
735	<	+ t) \right\rangle_{t_0} dt,
736	<	\end{equation}
737	<	although this method converges extremely slowly and is not practical
738	<	for obtaining viscosities from molecular dynamics simulations.
883	>	which converges much more rapidly in molecular dynamics simulations
884	>	than the traditional Green-Kubo formula.
885
886		The Langevin dynamics for the different model systems were performed
887		at the same temperature as the average temperature of the
#	Line 761 | Line 907 | have used in this study, there are differences between
907		compute the diffusive behavior for motion parallel to each body-fixed
908		axis by projecting the displacement of the particle onto the
909		body-fixed reference frame at $t=0$.  With an isotropic solvent, as we
910	<	have used in this study, there are differences between the three
911	<	diffusion constants, but these must converge to the same value at
912	<	longer times.  Translational diffusion constants for the different
913	<	shaped models are shown in table \ref{tab:translation}.
910	>	have used in this study, there may be differences between the three
911	>	diffusion constants at short times, but these must converge to the
912	>	same value at longer times.  Translational diffusion constants for the
913	>	different shaped models are shown in table \ref{tab:translation}.
914
915		In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational
916		diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object
#	Line 844 | Line 990 | hydrodynamic flows is well known, giving translational
990		an arbitrary value of 0.8 kcal/mol.
991
992		The Stokes-Einstein behavior of large spherical particles in
993	<	hydrodynamic flows is well known, giving translational friction
994	<	coefficients of $6 \pi \eta R$ (stick boundary conditions) and
995	<	rotational friction coefficients of $8 \pi \eta R^3$.  Recently,
996	<	Schmidt and Skinner have computed the behavior of spherical tag
997	<	particles in molecular dynamics simulations, and have shown that {\it
998	<	slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more
999	<	appropriate for molecule-sized spheres embedded in a sea of spherical
1000	<	solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx}
993	>	hydrodynamic flows with ``stick'' boundary conditions is well known,
994	>	and is given in Eqs. (\ref{eq:StokesTranslation}) and
995	>	(\ref{eq:StokesRotation}).  Recently, Schmidt and Skinner have
996	>	computed the behavior of spherical tag particles in molecular dynamics
997	>	simulations, and have shown that {\it slip} boundary conditions
998	>	($\Xi_{tt} = 4 \pi \eta \rho$) may be more appropriate for
999	>	molecule-sized spheres embedded in a sea of spherical solvent
1000	>	particles.\cite{Schmidt:2004fj,Schmidt:2003kx}
1001
1002		Our simulation results show similar behavior to the behavior observed
1003		by Schmidt and Skinner.  The diffusion constant obtained from our
#	Line 876 | Line 1022 | D = \frac{k_B T}{6 \pi \eta a} G(\rho),
1022		can be combined to give a single translational diffusion
1023		constant,\cite{Berne90}
1024		\begin{equation}
1025	<	D = \frac{k_B T}{6 \pi \eta a} G(\rho),
1025	>	D = \frac{k_B T}{6 \pi \eta a} G(s),
1026		\label{Dperrin}
1027		\end{equation}
1028		as well as a single rotational diffusion coefficient,
1029		\begin{equation}
1030	<	\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2)
1031	<	G(\rho) - 1}{1 - \rho^4} \right\}.
1030	>	\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - s^2)
1031	>	G(s) - 1}{1 - s^4} \right\}.
1032		\label{ThetaPerrin}
1033		\end{equation}
1034	<	In these expressions, $G(\rho)$ is a function of the axial ratio
1035	<	($\rho = b / a$), which for prolate ellipsoids, is
1034	>	In these expressions, $G(s)$ is a function of the axial ratio
1035	>	($s = b / a$), which for prolate ellipsoids, is
1036		\begin{equation}
1037	<	G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 -
892	<	\rho^2)^{1/2}}{\rho} \right\}
1037	>	G(s) = (1- s^2)^{-1/2} \ln \left\{ \frac{1 + (1 - s^2)^{1/2}}{s} \right\}
1038		\label{GPerrin}
1039		\end{equation}
1040		Again, there is some uncertainty about the correct boundary conditions
#	Line 914 | Line 1059 | The translational diffusion constants from the microca
1059		exact treatment of the diffusion tensor as well as the rough-shell
1060		model for the ellipsoid.
1061
1062	<	The translational diffusion constants from the microcanonical simulations
1063	<	agree well with the predictions of the Perrin model, although the rotational
1064	<	correlation times are a factor of 2 shorter than expected from hydrodynamic
1065	<	theory.  One explanation for the slower rotation
1066	<	of explicitly-solvated ellipsoids is the possibility that solute-solvent
1067	<	collisions happen at both ends of the solute whenever the principal
1068	<	axis of the ellipsoid is turning. In the upper portion of figure
1069	<	\ref{fig:explanation} we sketch a physical picture of this explanation.
1070	<	Since our Langevin integrator is providing nearly quantitative agreement with
1071	<	the Perrin model, it also predicts orientational diffusion for ellipsoids that
1072	<	exceed explicitly solvated correlation times by a factor of two.
1062	>	The translational diffusion constants from the microcanonical
1063	>	simulations agree well with the predictions of the Perrin model,
1064	>	although the {\it rotational} correlation times are a factor of 2
1065	>	shorter than expected from hydrodynamic theory.  One explanation for
1066	>	the slower rotation of explicitly-solvated ellipsoids is the
1067	>	possibility that solute-solvent collisions happen at both ends of the
1068	>	solute whenever the principal axis of the ellipsoid is turning. In the
1069	>	upper portion of figure \ref{fig:explanation} we sketch a physical
1070	>	picture of this explanation.  Since our Langevin integrator is
1071	>	providing nearly quantitative agreement with the Perrin model, it also
1072	>	predicts orientational diffusion for ellipsoids that exceed explicitly
1073	>	solvated correlation times by a factor of two.
1074
1075		\subsection{Rigid dumbbells}
1076		Perhaps the only {\it composite} rigid body for which analytic
#	Line 932 | Line 1078 | model. Equation (\ref{introEquation:oseenTensor}) abov
1078		two-sphere dumbbell model.  This model consists of two non-overlapping
1079		spheres held by a rigid bond connecting their centers. There are
1080		competing expressions for the 6x6 resistance tensor for this
1081	<	model. Equation (\ref{introEquation:oseenTensor}) above gives the
1082	<	original Oseen tensor, while the second order expression introduced by
1083	<	Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la
938	<	Torre and Bloomfield,\cite{Torre1977} is given above as
1081	>	model. The second order expression introduced by Rotne and
1082	>	Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la Torre and
1083	>	Bloomfield,\cite{Torre1977} is given above as
1084		Eq. (\ref{introEquation:RPTensorNonOverlapped}).  In our case, we use
1085		a model dumbbell in which the two spheres are identical Lennard-Jones
1086		particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at
#	Line 948 | Line 1093 | those derived from the 6 x 6 tensors mentioned above).
1093		motion in a flow {\it perpendicular} to the inter-sphere
1094		axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it
1095		orientational} correlation times for this model system (other than
1096	<	those derived from the 6 x 6 tensors mentioned above).
1096	>	those derived from the 6 x 6 tensor mentioned above).
1097
1098		The bead model for this model system comprises the two large spheres
1099		by themselves, while the rough shell approximation used 3368 separate
#	Line 1016 | Line 1161 | rotational diffusion.
1161		} \label{fig:explanation}
1162		\end{figure}
1163
1019	–
1020	–
1164		\subsection{Composite banana-shaped molecules}
1165		Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids have
1166		been used by Orlandi {\it et al.} to observe mesophases in
#	Line 1029 | Line 1172 | A reference system composed of a single banana rigid b
1172		behavior of this model, we have left out the dipolar interactions of
1173		the original Orlandi model.
1174
1175	<	A reference system composed of a single banana rigid body embedded in a
1176	<	sea of 1929 solvent particles was created and run under standard
1177	<	(microcanonical) molecular dynamics.  The resulting viscosity of this
1178	<	mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})).
1179	<	To calculate the hydrodynamic properties of the banana rigid body model,
1180	<	we created a rough shell (see Fig.~\ref{fig:roughShell}), in which
1181	<	the banana is represented as a ``shell'' made of 3321 identical beads
1182	<	(0.25 \AA\  in diameter) distributed on the surface.  Applying the
1183	<	procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
1184	<	identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 \AA), as
1185	<	well as the resistance tensor,
1186	<	\begin{equation*}
1044	<	\Xi =
1045	<	\left( {\begin{array}{*{20}c}
1046	<	0.9261 & 0 & 0&0&0.08585&0.2057\\
1047	<	0& 0.9270&-0.007063& 0.08585&0&0\\
1048	<	0&-0.007063&0.7494&0.2057&0&0\\
1049	<	0&0.0858&0.2057& 58.64& 0&0\\0.08585&0&0&0&48.30&3.219&\\0.2057&0&0&0&3.219&10.7373\\\end{array}} \right),
1050	<	\end{equation*}
1051	<	where the units for translational, translation-rotation coupling and
1052	<	rotational tensors are (kcal fs / mol \AA$^2$), (kcal fs / mol \AA\ rad),
1053	<	and (kcal fs / mol rad$^2$), respectively.
1175	>	A reference system composed of a single banana rigid body embedded in
1176	>	a sea of 1929 solvent particles was created and run under standard
1177	>	(microcanonical) molecular dynamics.  The resulting viscosity of this
1178	>	mixture was 0.298 centipoise (as estimated using
1179	>	Eq. (\ref{eq:shear})).  To calculate the hydrodynamic properties of
1180	>	the banana rigid body model, we created a rough shell (see
1181	>	Fig.~\ref{fig:roughShell}), in which the banana is represented as a
1182	>	``shell'' made of 3321 identical beads (0.25 \AA\ in diameter)
1183	>	distributed on the surface.  Applying the procedure described in
1184	>	Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
1185	>	identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0
1186	>	\AA).
1187
1188	<	The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor)
1189	<	are essentially quantitative for translational diffusion of this model.
1190	<	Orientational correlation times under the Langevin rigid-body integrator
1191	<	are within 11\% of the values obtained from explicit solvent, but these
1192	<	models also exhibit some solvent inaccessible surface area in the
1193	<	explicitly-solvated case.
1188	>	The Langevin rigid-body integrator (and the hydrodynamic diffusion
1189	>	tensor) are essentially quantitative for translational diffusion of
1190	>	this model.  Orientational correlation times under the Langevin
1191	>	rigid-body integrator are within 11\% of the values obtained from
1192	>	explicit solvent, but these models also exhibit some solvent
1193	>	inaccessible surface area in the explicitly-solvated case.
1194
1195		\subsection{Composite sphero-ellipsoids}
1063	–	Spherical heads perched on the ends of Gay-Berne ellipsoids have been
1064	–	used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01}
1196
1197	<	MORE DETAILS
1197	>	Spherical heads perched on the ends of Gay-Berne ellipsoids have been
1198	>	used recently as models for lipid
1199	>	molecules.\cite{SunX._jp0762020,Ayton01} A reference system composed
1200	>	of a single lipid rigid body embedded in a sea of 1929 solvent
1201	>	particles was created and run under a microcanonical ensemble.  The
1202	>	resulting viscosity of this mixture was 0.349 centipoise (as estimated
1203	>	using Eq. (\ref{eq:shear})).  To calculate the hydrodynamic properties
1204	>	of the lipid rigid body model, we created a rough shell (see
1205	>	Fig.~\ref{fig:roughShell}), in which the lipid is represented as a
1206	>	``shell'' made of 3550 identical beads (0.25 \AA\ in diameter)
1207	>	distributed on the surface.  Applying the procedure described by
1208	>	Eq. (\ref{introEquation:ResistanceTensorArbitraryOrigin}), we
1209	>	identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46
1210	>	\AA).
1211
1212	+	The translational diffusion constants and rotational correlation times
1213	+	obtained using the Langevin rigid-body integrator (and the
1214	+	hydrodynamic tensor) are essentially quantitative when compared with
1215	+	the explicit solvent simulations for this model system.
1216
1217	<	\subsection{Summary}
1218	<	According to our simulations, the langevin dynamics is a reliable
1219	<	theory to apply to replace the explicit solvents, especially for the
1220	<	translation properties. For large molecules, the rotation properties
1221	<	are also mimiced reasonablly well.
1217	>	\subsection{Summary of comparisons with explicit solvent simulations}
1218	>	The Langevin rigid-body integrator we have developed is a reliable way
1219	>	to replace explicit solvent simulations in cases where the detailed
1220	>	solute-solvent interactions do not greatly impact the behavior of the
1221	>	solute.  As such, it has the potential to greatly increase the length
1222	>	and time scales of coarse grained simulations of large solvated
1223	>	molecules.  In cases where the dielectric screening of the solvent, or
1224	>	specific solute-solvent interactions become important for structural
1225	>	or dynamic features of the solute molecule, this integrator may be
1226	>	less useful.  However, for the kinds of coarse-grained modeling that
1227	>	have become popular in recent years (ellipsoidal side chains, rigid
1228	>	bodies, and molecular-scale models), this integrator may prove itself
1229	>	to be quite valuable.
1230
1231	+	\begin{figure}
1232	+	\centering
1233	+	\includegraphics[width=\linewidth]{graph}
1234	+	\caption[Mean squared displacements and orientational
1235	+	correlation functions for each of the model rigid bodies.]{The
1236	+	mean-squared displacements ($\langle r^2(t) \rangle$) and
1237	+	orientational correlation functions ($C_2(t)$) for each of the model
1238	+	rigid bodies studied.  The circles are the results for microcanonical
1239	+	simulations with explicit solvent molecules, while the other data sets
1240	+	are results for Langevin dynamics using the different hydrodynamic
1241	+	tensor approximations.  The Perrin model for the ellipsoids is
1242	+	considered the ``exact'' hydrodynamic behavior (this can also be said
1243	+	for the translational motion of the dumbbell operating under the bead
1244	+	model). In most cases, the various hydrodynamics models reproduce
1245	+	each other quantitatively.}
1246	+	\label{fig:results}
1247	+	\end{figure}
1248	+
1249		\begin{table*}
1250		\begin{minipage}{\linewidth}
1251		\begin{center}
1252		\caption{Translational diffusion constants (D) for the model systems
1253		calculated using microcanonical simulations (with explicit solvent),
1254		theoretical predictions, and Langevin simulations (with implicit solvent).
1255	<	Analytical solutions for the exactly-solved hydrodynamics models are
1256	<	from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936}
1255	>	Analytical solutions for the exactly-solved hydrodynamics models are obtained
1256	>	from: Stokes' law (sphere), and Refs. \citen{Perrin1934} and \citen{Perrin1936}
1257		(ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq}
1258		(dumbbell). The other model systems have no known analytic solution.
1259	<	All  diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (=
1259	>	All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (=
1260		$10^{-4}$ \AA$^2$  / fs). }
1261		\begin{tabular}{lccccccc}
1262		\hline
#	Line 1090 | Line 1264 | sphere    & 0.261  & ?    & & 2.59 & exact       & 2.5
1264		\cline{2-3} \cline{5-7}
1265		model & $\eta$ (centipoise)  & D & & Analytical & method & Hydrodynamics & simulation \\
1266		\hline
1267	<	sphere    & 0.261  & ?    & & 2.59 & exact       & 2.59 & 2.56 \\
1267	>	sphere    & 0.279  & 3.06 & & 2.42 & exact       & 2.42 & 2.33 \\
1268		ellipsoid & 0.255  & 2.44 & & 2.34 & exact       & 2.34 & 2.37 \\
1269		& 0.255  & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\
1270	<	dumbbell  & 0.322  & ?    & & 1.57 & bead model  & 1.57 & 1.57 \\
1271	<	& 0.322  & ?    & & 1.57 & rough shell & ?    & ?    \\
1270	>	dumbbell  & 0.308  & 2.06 & & 1.64 & bead model  & 1.65 & 1.62 \\
1271	>	& 0.308  & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\
1272		banana    & 0.298  & 1.53 & &      & rough shell & 1.56 & 1.55 \\
1273	<	lipid     & 0.349  & 0.96 & &      & rough shell & 1.33 & 1.32 \\
1273	>	lipid     & 0.349  & 1.41 & &      & rough shell & 1.33 & 1.32 \\
1274		\end{tabular}
1275		\label{tab:translation}
1276		\end{center}
#	Line 1119 | Line 1293 | sphere    & 0.261  &      & & 9.06 & exact       & 9.0
1293		\cline{2-3} \cline{5-7}
1294		model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic  & simulation \\
1295		\hline
1296	<	sphere    & 0.261  &      & & 9.06 & exact       & 9.06 & 9.11 \\
1296	>	sphere    & 0.279  &      & & 9.69 & exact       & 9.69 & 9.64 \\
1297		ellipsoid & 0.255  & 46.7 & & 22.0 & exact       & 22.0 & 22.2 \\
1298		& 0.255  & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\
1299	<	dumbbell  & 0.322  & 14.0 & &      & bead model  & 52.3 & 52.8 \\
1300	<	& 0.322  & 14.0 & &      & rough shell & ?    & ?    \\
1299	>	dumbbell  & 0.308  & 14.1 & &      & bead model  & 50.0 & 50.1 \\
1300	>	& 0.308  & 14.1 & &      & rough shell & 41.5 & 41.3 \\
1301		banana    & 0.298  & 63.8 & &      & rough shell & 70.9 & 70.9 \\
1302		lipid     & 0.349  & 78.0 & &      & rough shell & 76.9 & 77.9 \\
1303		\hline
#	Line 1135 | Line 1309 | The Langevin dynamics integrator was applied to study
1309
1310		\section{Application: A rigid-body lipid bilayer}
1311
1312	<	The Langevin dynamics integrator was applied to study the formation of
1313	<	corrugated structures emerging from simulations of the coarse grained
1314	<	lipid molecular models presented above.  The initial configuration is
1315	<	taken from our molecular dynamics studies on lipid bilayers with
1316	<	lennard-Jones sphere solvents. The solvent molecules were excluded
1317	<	from the system, and the experimental value for the viscosity of water
1318	<	at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects
1319	<	of the solvent.  The absence of explicit solvent molecules and the
1320	<	stability of the integrator allowed us to take timesteps of 50 fs.  A
1321	<	total simulation run time of 100 ns was sampled.
1322	<	Fig. \ref{fig:bilayer} shows the configuration of the system after 100
1323	<	ns, and the ripple structure remains stable during the entire
1324	<	trajectory.  Compared with using explicit bead-model solvent
1325	<	molecules, the efficiency of the simulation has increased by an order
1326	<	of magnitude.
1312	>	To test the accuracy and efficiency of the new integrator, we applied
1313	>	it to study the formation of corrugated structures emerging from
1314	>	simulations of the coarse grained lipid molecular models presented
1315	>	above.  The initial configuration is taken from earlier molecular
1316	>	dynamics studies on lipid bilayers which had used spherical
1317	>	(Lennard-Jones) solvent particles and moderate (480 solvated lipid
1318	>	molecules) system sizes.\cite{SunX._jp0762020} the solvent molecules
1319	>	were excluded from the system and the box was replicated three times
1320	>	in the x- and y- axes (giving a single simulation cell containing 4320
1321	>	lipids).  The experimental value for the viscosity of water at 20C
1322	>	($\eta = 1.00$ cp) was used with the Langevin integrator to mimic the
1323	>	hydrodynamic effects of the solvent.  The absence of explicit solvent
1324	>	molecules and the stability of the integrator allowed us to take
1325	>	timesteps of 50 fs. A simulation run time of 30 ns was sampled to
1326	>	calculate structural properties.  Fig. \ref{fig:bilayer} shows the
1327	>	configuration of the system after 30 ns.  Structural properties of the
1328	>	bilayer (e.g. the head and body $P_2$ order parameters) are nearly
1329	>	identical to those obtained via solvated molecular dynamics. The
1330	>	ripple structure remained stable during the entire trajectory.
1331	>	Compared with using explicit bead-model solvent molecules, the 30 ns
1332	>	trajectory for 4320 lipids with the Langevin integrator is now {\it
1333	>	faster} on the same hardware than the same length trajectory was for
1334	>	the 480-lipid system previously studied.
1335
1336		\begin{figure}
1337		\centering
1338		\includegraphics[width=\linewidth]{bilayer}
1339		\caption[Snapshot of a bilayer of rigid-body models for lipids]{A
1340	<	snapshot of a bilayer composed of rigid-body models for lipid
1340	>	snapshot of a bilayer composed of 4320 rigid-body models for lipid
1341		molecules evolving using the Langevin integrator described in this
1342		work.} \label{fig:bilayer}
1343		\end{figure}
1344
1345		\section{Conclusions}
1346
1347	<	We have presented a new Langevin algorithm by incorporating the
1348	<	hydrodynamics properties of arbitrary shaped molecules into an
1349	<	advanced symplectic integration scheme. Further studies in systems
1350	<	involving banana shaped molecules illustrated that the dynamic
1351	<	properties could be preserved by using this new algorithm as an
1352	<	implicit solvent model.
1347	>	We have presented a new algorithm for carrying out Langevin dynamics
1348	>	simulations on complex rigid bodies by incorporating the hydrodynamic
1349	>	resistance tensors for arbitrary shapes into a stable and efficient
1350	>	integration scheme.  The integrator gives quantitative agreement with
1351	>	both analytic and approximate hydrodynamic theories, and works
1352	>	reasonably well at reproducing the solute dynamical properties
1353	>	(diffusion constants, and orientational relaxation times) from
1354	>	explicitly-solvated simulations.  For the cases where there are
1355	>	discrepancies between our Langevin integrator and the explicit solvent
1356	>	simulations, two features of molecular simulations help explain the
1357	>	differences.
1358
1359	+	First, the use of ``stick'' boundary conditions for molecular-sized
1360	+	solutes in a sea of similarly-sized solvent particles may be
1361	+	problematic.  We are certainly not the first group to notice this
1362	+	difference between hydrodynamic theories and explicitly-solvated
1363	+	molecular
1364	+	simulations.\cite{Schmidt:2004fj,Schmidt:2003kx,Ravichandran:1999fk,TANG:1993lr}
1365	+	The problem becomes particularly noticable in both the translational
1366	+	diffusion of the spherical particles and the rotational diffusion of
1367	+	the ellipsoids.  In both of these cases it is clear that the
1368	+	approximations that go into hydrodynamics are the source of the error,
1369	+	and not the integrator itself.
1370
1371	+	Second, in the case of structures which have substantial surface area
1372	+	that is inaccessible to solvent particles, the hydrodynamic theories
1373	+	(and the Langevin integrator) may overestimate the effects of solvent
1374	+	friction because they overestimate the exposed surface area of the
1375	+	rigid body.  This is particularly noticable in the rotational
1376	+	diffusion of the dumbbell model.  We believe that given a solvent of
1377	+	known radius, it may be possible to modify the rough shell approach to
1378	+	place beads on solvent-accessible surface, instead of on the geometric
1379	+	surface defined by the van der Waals radii of the components of the
1380	+	rigid body.  Further work to confirm the behavior of this new
1381	+	approximation is ongoing.
1382	+
1383		\section{Acknowledgments}
1384		Support for this project was provided by the National Science
1385		Foundation under grant CHE-0134881. T.L. also acknowledges the
1386	<	financial support from center of applied mathematics at University
1387	<	of Notre Dame.
1386	>	financial support from Center of Applied Mathematics at University of
1387	>	Notre Dame.
1388	>
1389	>	\end{doublespace}
1390		\newpage
1391
1392	<	\bibliographystyle{jcp}
1392	>	\bibliographystyle{jcp2}
1393		\bibliography{langevin}
1182	–
1394		\end{document}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines