ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/langevin/langevin.tex
Revision: 3340
Committed: Thu Jan 31 22:13:55 2008 UTC (16 years, 7 months ago) by gezelter
Content type: application/x-tex
File size: 58317 byte(s)
Log Message:
edits

File Contents

# User Rev Content
1 tim 2746 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2     %\documentclass[aps,prb,preprint]{revtex4}
3     \documentclass[11pt]{article}
4     \usepackage{endfloat}
5 gezelter 3337 \usepackage{amsmath}
6 tim 2746 \usepackage{amssymb}
7     \usepackage{times}
8     \usepackage{mathptmx}
9     \usepackage{setspace}
10     \usepackage{tabularx}
11     \usepackage{graphicx}
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}
20 gezelter 3299 \renewcommand\citemid{\ } % no comma in optional referenc note
21 tim 2746
22     \begin{document}
23    
24 gezelter 3205 \title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape }
25 tim 2746
26 gezelter 3299 \author{Xiuquan Sun, Teng Lin and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
27 tim 2746 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
35    
36     \begin{abstract}
37    
38     \end{abstract}
39    
40     \newpage
41    
42     %\narrowtext
43    
44     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45     % BODY OF TEXT
46     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47    
48     \section{Introduction}
49    
50     %applications of langevin dynamics
51 gezelter 3316 Langevin dynamics, which mimics a simple heat bath with stochastic and
52     dissipative forces, has been applied in a variety of situations as an
53     alternative to molecular dynamics with explicit solvent molecules.
54     The stochastic treatment of the solvent allows the use of simulations
55     with substantially longer time and length scales. In general, the
56     dynamic and structural properties obtained from Langevin simulations
57     agree quite well with similar properties obtained from explicit
58     solvent simulations.
59    
60     Recent examples of the usefulness of Langevin simulations include a
61     study of met-enkephalin in which Langevin simulations predicted
62     dynamical properties that were largely in agreement with explicit
63     solvent simulations.\cite{Shen2002} By applying Langevin dynamics with
64     the UNRES model, Liow and his coworkers suggest that protein folding
65     pathways can be explored within a reasonable amount of
66     time.\cite{Liwo2005}
67    
68     The stochastic nature of Langevin dynamics also enhances the sampling
69     of the system and increases the probability of crossing energy
70     barriers.\cite{Cui2003,Banerjee2004} Combining Langevin dynamics with
71 gezelter 3333 Kramers' theory, Klimov and Thirumalai identified free-energy
72 gezelter 3316 barriers by studying the viscosity dependence of the protein folding
73     rates.\cite{Klimov1997} In order to account for solvent induced
74     interactions missing from the implicit solvent model, Kaya
75     incorporated a desolvation free energy barrier into protein
76     folding/unfolding studies and discovered a higher free energy barrier
77 xsun 3317 between the native and denatured states.\cite{HuseyinKaya07012005}
78 gezelter 3316
79     Because of its stability against noise, Langevin dynamics has also
80     proven useful for studying remagnetization processes in various
81     systems.\cite{Palacios1998,Berkov2002,Denisov2003} [Check: For
82 tim 2746 instance, the oscillation power spectrum of nanoparticles from
83 gezelter 3316 Langevin dynamics has the same peak frequencies for different wave
84     vectors, which recovers the property of magnetic excitations in small
85     finite structures.\cite{Berkov2005a}]
86 tim 2746
87 gezelter 3316 In typical LD simulations, the friction and random forces on
88 gezelter 3333 individual atoms are taken from Stokes' law,
89 gezelter 3316 \begin{eqnarray}
90     m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\
91     \langle R(t) \rangle & = & 0 \\
92     \langle R(t) R(t') \rangle & = & 2 k_B T \xi m \delta(t - t')
93     \end{eqnarray}
94     where $\xi \approx 6 \pi \eta a$. Here $\eta$ is the viscosity of the
95     implicit solvent, and $a$ is the hydrodynamic radius of the atom.
96 tim 2746
97 gezelter 3333 The use of rigid substructures,\cite{Chun:2000fj}
98     coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunGezelter08}
99     and ellipsoidal representations of protein side chains~\cite{Fogolari:1996lr}
100     has made the use of the Stokes-Einstein approximation problematic. A
101     rigid substructure moves as a single unit with orientational as well
102     as translational degrees of freedom. This requires a more general
103 gezelter 3316 treatment of the hydrodynamics than the spherical approximation
104     provides. The atoms involved in a rigid or coarse-grained structure
105     should properly have solvent-mediated interactions with each
106     other. The theory of interactions {\it between} bodies moving through
107     a fluid has been developed over the past century and has been applied
108     to simulations of Brownian
109 gezelter 3333 motion.\cite{FIXMAN:1986lr,Ramachandran1996}
110 tim 2746
111 gezelter 3333 In order to account for the diffusion anisotropy of arbitrarily-shaped
112     particles, Fernandes and Garc\'{i}a de la Torre improved the original
113     Brownian dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by
114     incorporating a generalized $6\times6$ diffusion tensor and
115     introducing a rotational evolution scheme consisting of three
116     consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are
117     introduced into the system due to the arbitrary order of applying the
118     noncommuting rotation operators.\cite{Beard2003} Based on the
119     observation the momentum relaxation time is much less than the time
120     step, one may ignore the inertia in Brownian dynamics. However, the
121     assumption of zero average acceleration is not always true for
122     cooperative motion which is common in proteins. An inertial Brownian
123     dynamics (IBD) was proposed to address this issue by adding an
124     inertial correction term.\cite{Beard2000} As a complement to IBD which
125     has a lower bound in time step because of the inertial relaxation
126     time, long-time-step inertial dynamics (LTID) can be used to
127     investigate the inertial behavior of linked polymer segments in a low
128     friction regime.\cite{Beard2000} LTID can also deal with the
129     rotational dynamics for nonskew bodies without translation-rotation
130     coupling by separating the translation and rotation motion and taking
131     advantage of the analytical solution of hydrodynamics
132     properties. However, typical nonskew bodies like cylinders and
133     ellipsoids are inadequate to represent most complex macromolecular
134     assemblies. There is therefore a need for incorporating the
135     hydrodynamics of complex (and potentially skew) rigid bodies in the
136     library of methods available for performing Langevin simulations.
137    
138 gezelter 3316 \subsection{Rigid Body Dynamics}
139     Rigid bodies are frequently involved in the modeling of large
140     collections of particles that move as a single unit. In molecular
141     simulations, rigid bodies have been used to simplify protein-protein
142 gezelter 3333 docking,\cite{Gray2003} and lipid bilayer
143     simulations.\cite{SunGezelter08} Many of the water models in common
144     use are also rigid-body
145     models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are
146     typically evolved using constraints rather than rigid body equations
147     of motion.
148 gezelter 3316
149 gezelter 3333 Euler angles are a natural choice to describe the rotational degrees
150 gezelter 3337 of freedom. However, due to $\frac{1}{\sin \theta}$ singularities, the
151 gezelter 3333 numerical integration of corresponding equations of these motion can
152     become inaccurate (and inefficient). Although the use of multiple
153     sets of Euler angles can overcome this problem,\cite{Barojas1973} the
154     computational penalty and the loss of angular momentum conservation
155     remain. A singularity-free representation utilizing quaternions was
156     developed by Evans in 1977.\cite{Evans1977} The Evans quaternion
157     approach uses a nonseparable Hamiltonian, and this has prevented
158     symplectic algorithms from being utilized until very
159     recently.\cite{Miller2002}
160 gezelter 3316
161 gezelter 3333 Another approach is the application of holonomic constraints to the
162     atoms belonging to the rigid body. Each atom moves independently
163     under the normal forces deriving from potential energy and constraints
164     are used to guarantee rigidity. However, due to their iterative
165     nature, the SHAKE and RATTLE algorithms converge very slowly when the
166     number of constraints (and the number of particles that belong to the
167     rigid body) increases.\cite{Ryckaert1977,Andersen1983}
168 tim 2746
169 gezelter 3333 In order to develop a stable and efficient integration scheme that
170     preserves most constants of the motion, symplectic propagators are
171     necessary. By introducing a conjugate momentum to the rotation matrix
172     $Q$ and re-formulating Hamilton's equations, a symplectic
173     orientational integrator, RSHAKE,\cite{Kol1997} was proposed to evolve
174     rigid bodies on a constraint manifold by iteratively satisfying the
175     orthogonality constraint $Q^T Q = 1$. An alternative method using the
176     quaternion representation was developed by Omelyan.\cite{Omelyan1998}
177     However, both of these methods are iterative and suffer from some
178     related inefficiencies. A symplectic Lie-Poisson integrator for rigid
179     bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes
180     most of the limitations mentioned above and is therefore the basis for
181     our Langevin integrator.
182 gezelter 3316
183 tim 2746 The goal of the present work is to develop a Langevin dynamics
184 tim 2999 algorithm for arbitrary-shaped rigid particles by integrating the
185 gezelter 3316 accurate estimation of friction tensor from hydrodynamics theory into
186     a symplectic rigid body dynamics propagator. In the sections below,
187 gezelter 3333 we review some of the theory of hydrodynamic tensors developed
188     primarily for Brownian simulations of multi-particle systems, we then
189     present our integration method for a set of generalized Langevin
190     equations of motion, and we compare the behavior of the new Langevin
191     integrator to dynamical quantities obtained via explicit solvent
192     molecular dynamics.
193 tim 2746
194 gezelter 3316 \subsection{\label{introSection:frictionTensor}The Friction Tensor}
195     Theoretically, a complete friction kernel can be determined using the
196 tim 2999 velocity autocorrelation function. However, this approach becomes
197 gezelter 3333 impractical when the solute becomes complex. Instead, various
198 gezelter 3316 approaches based on hydrodynamics have been developed to calculate the
199     friction coefficients. In general, the friction tensor $\Xi$ is a
200     $6\times 6$ matrix given by
201     \begin{equation}
202 gezelter 3333 \Xi = \left( \begin{array}{*{20}c}
203     \Xi^{tt} & \Xi^{rt} \\
204     \Xi^{tr} & \Xi^{rr} \\
205     \end{array} \right).
206 gezelter 3316 \end{equation}
207     Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and
208     rotational resistance (friction) tensors respectively, while
209     $\Xi^{tr}$ is translation-rotation coupling tensor and $\Xi^{rt}$ is
210     rotation-translation coupling tensor. When a particle moves in a
211 xsun 3338 fluid, it may experience friction force ($\mathbf{f}_f$) and torque
212 gezelter 3316 ($\mathbf{\tau}_f$) in opposition to the directions of the velocity
213     ($\mathbf{v}$) and body-fixed angular velocity ($\mathbf{\omega}$),
214     \begin{equation}
215 tim 2746 \left( \begin{array}{l}
216 xsun 3338 \mathbf{f}_f \\
217 gezelter 3316 \mathbf{\tau}_f \\
218 gezelter 3333 \end{array} \right) = - \left( \begin{array}{*{20}c}
219     \Xi^{tt} & \Xi^{rt} \\
220     \Xi^{tr} & \Xi^{rr} \\
221     \end{array} \right)\left( \begin{array}{l}
222 gezelter 3316 \mathbf{v} \\
223     \mathbf{\omega} \\
224     \end{array} \right).
225     \end{equation}
226 tim 2746
227 tim 2999 \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
228 gezelter 3316 For a spherical particle under ``stick'' boundary conditions, the
229     translational and rotational friction tensors can be calculated from
230 gezelter 3333 Stokes' law,
231 gezelter 3316 \begin{equation}
232 gezelter 3333 \Xi^{tt} = \left( \begin{array}{*{20}c}
233 tim 2999 {6\pi \eta R} & 0 & 0 \\
234     0 & {6\pi \eta R} & 0 \\
235     0 & 0 & {6\pi \eta R} \\
236 gezelter 3333 \end{array} \right)
237 gezelter 3316 \end{equation}
238 tim 2999 and
239 gezelter 3316 \begin{equation}
240 gezelter 3333 \Xi^{rr} = \left( \begin{array}{*{20}c}
241 tim 2999 {8\pi \eta R^3 } & 0 & 0 \\
242     0 & {8\pi \eta R^3 } & 0 \\
243     0 & 0 & {8\pi \eta R^3 } \\
244 gezelter 3333 \end{array} \right)
245 gezelter 3316 \end{equation}
246 tim 2999 where $\eta$ is the viscosity of the solvent and $R$ is the
247     hydrodynamic radius.
248    
249     Other non-spherical shapes, such as cylinders and ellipsoids, are
250 gezelter 3316 widely used as references for developing new hydrodynamics theories,
251 tim 2999 because their properties can be calculated exactly. In 1936, Perrin
252 gezelter 3333 extended Stokes' law to general ellipsoids which are given in
253     Cartesian coordinates by~\cite{Perrin1934,Perrin1936}
254 gezelter 3316 \begin{equation}
255 gezelter 3333 \frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1.
256 gezelter 3316 \end{equation}
257 gezelter 3333 Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the
258     complexity of the elliptic integral, only uniaxial ellipsoids, either
259     prolate ($a \ge b = c$) or oblate ($a < b = c$), can be solved
260     exactly. Introducing an elliptic integral parameter $S$ for prolate,
261 gezelter 3316 \begin{equation}
262     S = \frac{2}{\sqrt{a^2 - b^2}} \ln \frac{a + \sqrt{a^2 - b^2}}{b},
263     \end{equation}
264 gezelter 3333 and oblate,
265 gezelter 3316 \begin{equation}
266     S = \frac{2}{\sqrt {b^2 - a^2 }} \arctan \frac{\sqrt {b^2 - a^2}}{a},
267     \end{equation}
268 gezelter 3333 ellipsoids, one can write down the translational and rotational
269     resistance tensors:
270 tim 2999 \begin{eqnarray*}
271 gezelter 3316 \Xi_a^{tt} & = & 16\pi \eta \frac{a^2 - b^2}{(2a^2 - b^2 )S - 2a}. \\
272     \Xi_b^{tt} = \Xi_c^{tt} & = & 32\pi \eta \frac{a^2 - b^2 }{(2a^2 - 3b^2 )S + 2a},
273 tim 2999 \end{eqnarray*}
274 gezelter 3333 for oblate, and
275 tim 2999 \begin{eqnarray*}
276 gezelter 3316 \Xi_a^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^2 - b^2 )b^2}{2a - b^2 S}, \\
277     \Xi_b^{rr} = \Xi_c^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^4 - b^4)}{(2a^2 - b^2 )S - 2a}
278 tim 2999 \end{eqnarray*}
279 gezelter 3333 for prolate ellipsoids. For both spherical and ellipsoidal particles,
280     the translation-rotation and rotation-translation coupling tensors are
281 gezelter 3316 zero.
282 tim 2746
283 tim 2999 \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
284     Unlike spherical and other simply shaped molecules, there is no
285     analytical solution for the friction tensor for arbitrarily shaped
286     rigid molecules. The ellipsoid of revolution model and general
287     triaxial ellipsoid model have been used to approximate the
288 gezelter 3316 hydrodynamic properties of rigid bodies. However, the mapping from all
289     possible ellipsoidal spaces, $r$-space, to all possible combination of
290     rotational diffusion coefficients, $D$-space, is not
291     unique.\cite{Wegener1979} Additionally, because there is intrinsic
292     coupling between translational and rotational motion of rigid bodies,
293     general ellipsoids are not always suitable for modeling arbitrarily
294     shaped rigid molecules. A number of studies have been devoted to
295 tim 2999 determining the friction tensor for irregularly shaped rigid bodies
296 gezelter 3316 using more advanced methods where the molecule of interest was modeled
297     by a combinations of spheres\cite{Carrasco1999} and the hydrodynamics
298     properties of the molecule can be calculated using the hydrodynamic
299 gezelter 3333 interaction tensor.
300    
301     Consider a rigid assembly of $N$ beads immersed in a continuous
302     medium. Due to hydrodynamic interaction, the ``net'' velocity of $i$th
303     bead, $v'_i$ is different than its unperturbed velocity $v_i$,
304     \begin{equation}
305 tim 2999 v'_i = v_i - \sum\limits_{j \ne i} {T_{ij} F_j }
306 gezelter 3333 \end{equation}
307     where $F_i$ is the frictional force, and $T_{ij}$ is the hydrodynamic
308     interaction tensor. The frictional force on the $i^\mathrm{th}$ bead
309     is proportional to its ``net'' velocity
310 tim 2746 \begin{equation}
311 tim 2999 F_i = \zeta _i v_i - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
312     \label{introEquation:tensorExpression}
313 tim 2746 \end{equation}
314 tim 2999 This equation is the basis for deriving the hydrodynamic tensor. In
315     1930, Oseen and Burgers gave a simple solution to
316     Eq.~\ref{introEquation:tensorExpression}
317 tim 2746 \begin{equation}
318 tim 2999 T_{ij} = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
319     R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
320 tim 2746 \end{equation}
321 tim 2999 Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
322     A second order expression for element of different size was
323     introduced by Rotne and Prager\cite{Rotne1969} and improved by
324     Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977}
325 tim 2746 \begin{equation}
326 tim 2999 T_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
327     \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
328     _i^2 + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
329     \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
330     \label{introEquation:RPTensorNonOverlapped}
331 tim 2746 \end{equation}
332 tim 2999 Both of the Eq.~\ref{introEquation:oseenTensor} and
333     Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption
334     $R_{ij} \ge \sigma _i + \sigma _j$. An alternative expression for
335     overlapping beads with the same radius, $\sigma$, is given by
336 tim 2746 \begin{equation}
337 tim 2999 T_{ij} = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
338     \frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
339     \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
340     \label{introEquation:RPTensorOverlapped}
341 tim 2746 \end{equation}
342 tim 2999 To calculate the resistance tensor at an arbitrary origin $O$, we
343     construct a $3N \times 3N$ matrix consisting of $N \times N$
344     $B_{ij}$ blocks
345     \begin{equation}
346 gezelter 3333 B = \left( \begin{array}{*{20}c}
347     B_{11} & \ldots & B_{1N} \\
348 tim 2999 \vdots & \ddots & \vdots \\
349 gezelter 3333 B_{N1} & \cdots & B_{NN} \\
350     \end{array} \right),
351 tim 2999 \end{equation}
352     where $B_{ij}$ is given by
353 gezelter 3333 \begin{equation}
354 tim 2999 B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
355     )T_{ij}
356 gezelter 3333 \end{equation}
357 tim 2999 where $\delta _{ij}$ is the Kronecker delta function. Inverting the
358     $B$ matrix, we obtain
359 tim 2746 \[
360 tim 2999 C = B^{ - 1} = \left( {\begin{array}{*{20}c}
361     {C_{11} } & \ldots & {C_{1N} } \\
362     \vdots & \ddots & \vdots \\
363     {C_{N1} } & \cdots & {C_{NN} } \\
364     \end{array}} \right),
365 tim 2746 \]
366 tim 2999 which can be partitioned into $N \times N$ $3 \times 3$ block
367     $C_{ij}$. With the help of $C_{ij}$ and the skew matrix $U_i$
368 tim 2746 \[
369 tim 2999 U_i = \left( {\begin{array}{*{20}c}
370     0 & { - z_i } & {y_i } \\
371     {z_i } & 0 & { - x_i } \\
372     { - y_i } & {x_i } & 0 \\
373     \end{array}} \right)
374 tim 2746 \]
375 tim 2999 where $x_i$, $y_i$, $z_i$ are the components of the vector joining
376     bead $i$ and origin $O$, the elements of resistance tensor at
377     arbitrary origin $O$ can be written as
378     \begin{eqnarray}
379 xsun 3339 \label{introEquation:ResistanceTensorArbitraryOrigin}
380 tim 2999 \Xi _{}^{tt} & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
381     \Xi _{}^{tr} & = & \Xi _{}^{rt} = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
382 gezelter 3310 \Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } }
383     U_j + 6 \eta V {\bf I}. \notag
384 tim 2999 \end{eqnarray}
385 gezelter 3310 The final term in the expression for $\Xi^{rr}$ is correction that
386     accounts for errors in the rotational motion of certain kinds of bead
387     models. The additive correction uses the solvent viscosity ($\eta$)
388     as well as the total volume of the beads that contribute to the
389     hydrodynamic model,
390     \begin{equation}
391     V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3,
392     \end{equation}
393     where $\sigma_i$ is the radius of bead $i$. This correction term was
394     rigorously tested and compared with the analytical results for
395     two-sphere and ellipsoidal systems by Garcia de la Torre and
396     Rodes.\cite{Torre:1983lr}
397    
398    
399 tim 2999 The resistance tensor depends on the origin to which they refer. The
400     proper location for applying the friction force is the center of
401     resistance (or center of reaction), at which the trace of rotational
402     resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
403     Mathematically, the center of resistance is defined as an unique
404     point of the rigid body at which the translation-rotation coupling
405     tensors are symmetric,
406     \begin{equation}
407     \Xi^{tr} = \left( {\Xi^{tr} } \right)^T
408     \label{introEquation:definitionCR}
409     \end{equation}
410     From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
411     we can easily derive that the translational resistance tensor is
412     origin independent, while the rotational resistance tensor and
413     translation-rotation coupling resistance tensor depend on the
414     origin. Given the resistance tensor at an arbitrary origin $O$, and
415     a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
416     obtain the resistance tensor at $P$ by
417     \begin{equation}
418     \begin{array}{l}
419     \Xi _P^{tt} = \Xi _O^{tt} \\
420     \Xi _P^{tr} = \Xi _P^{rt} = \Xi _O^{tr} - U_{OP} \Xi _O^{tt} \\
421     \Xi _P^{rr} = \Xi _O^{rr} - U_{OP} \Xi _O^{tt} U_{OP} + \Xi _O^{tr} U_{OP} - U_{OP} \Xi _O^{{tr} ^{^T }} \\
422     \end{array}
423     \label{introEquation:resistanceTensorTransformation}
424     \end{equation}
425     where
426 tim 2746 \[
427 tim 2999 U_{OP} = \left( {\begin{array}{*{20}c}
428     0 & { - z_{OP} } & {y_{OP} } \\
429 xsun 3338 {z_{OP} } & 0 & { - x_{OP} } \\
430 tim 2999 { - y_{OP} } & {x_{OP} } & 0 \\
431     \end{array}} \right)
432 tim 2746 \]
433 tim 2999 Using Eq.~\ref{introEquation:definitionCR} and
434     Eq.~\ref{introEquation:resistanceTensorTransformation}, one can
435     locate the position of center of resistance,
436     \begin{eqnarray*}
437     \left( \begin{array}{l}
438     x_{OR} \\
439     y_{OR} \\
440     z_{OR} \\
441 gezelter 3333 \end{array} \right) & = &\left( \begin{array}{*{20}c}
442 tim 2999 {(\Xi _O^{rr} )_{yy} + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} } \\
443     { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz} + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} } \\
444     { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx} + (\Xi _O^{rr} )_{yy} } \\
445 gezelter 3333 \end{array} \right)^{ - 1} \\
446 tim 2999 & & \left( \begin{array}{l}
447     (\Xi _O^{tr} )_{yz} - (\Xi _O^{tr} )_{zy} \\
448     (\Xi _O^{tr} )_{zx} - (\Xi _O^{tr} )_{xz} \\
449     (\Xi _O^{tr} )_{xy} - (\Xi _O^{tr} )_{yx} \\
450     \end{array} \right) \\
451     \end{eqnarray*}
452 xsun 3338 where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector
453 tim 2999 joining center of resistance $R$ and origin $O$.
454 tim 2746
455    
456 gezelter 3310 \section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
457 gezelter 3337
458 tim 2999 Consider the Langevin equations of motion in generalized coordinates
459 tim 2746 \begin{equation}
460 gezelter 3337 \mathbf{M} \dot{\mathbf{V}}(t) = \mathbf{F}_{s}(t) +
461 xsun 3338 \mathbf{F}_{f}(t) + \mathbf{F}_{r}(t)
462 tim 2746 \label{LDGeneralizedForm}
463     \end{equation}
464 gezelter 3337 where $\mathbf{M}$ is a $6 \times 6$ diagonal mass matrix (which
465     includes the mass of the rigid body as well as the moments of inertia
466     in the body-fixed frame) and $\mathbf{V}$ is a generalized velocity,
467     $\mathbf{V} =
468     \left\{\mathbf{v},\mathbf{\omega}\right\}$. The right side of
469 gezelter 3333 Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a
470 gezelter 3337 system force $\mathbf{F}_{s}$, a frictional or dissipative force
471 xsun 3338 $\mathbf{F}_{f}$ and stochastic force $\mathbf{F}_{r}$. While the
472 gezelter 3340 evolution of the system in Newtonian mechanics is typically done in
473     the lab-fixed frame, it is convenient to handle the dynamics of rigid
474     bodies in the body-fixed frame. Thus the friction and random forces
475     are calculated in body-fixed frame and may be converted back to
476     lab-fixed frame using the rigid body's rotation matrix ($Q$):
477 gezelter 3333 \begin{equation}
478 gezelter 3340 \mathbf{F}_{f,r} =
479     \left( \begin{array}{c}
480     \mathbf{f}_{f,r} \\
481     \mathbf{\tau}_{f,r}
482     \end{array} \right)
483     =
484     \left( \begin{array}{c}
485     Q^{T} \mathbf{f}^{b}_{f,r} \\
486     Q^{T} \mathbf{\tau}^{b}_{f,r}
487     \end{array} \right)
488 gezelter 3333 \end{equation}
489 gezelter 3340 The body-fixed friction force, $\mathbf{F}_{f}^b$, is proportional to
490     the velocity at the center of resistance $\mathbf{v}_{R}^b$ (in the
491     body-fixed frame) and the angular velocity $\mathbf{\omega}$
492 tim 2746 \begin{equation}
493 gezelter 3337 \mathbf{F}_{f}^b (t) = \left( \begin{array}{l}
494     \mathbf{f}_{f}^b (t) \\
495     \mathbf{\tau}_{f}^b (t) \\
496 gezelter 3333 \end{array} \right) = - \left( \begin{array}{*{20}c}
497     \Xi_{R,t} & \Xi_{R,c}^T \\
498     \Xi_{R,c} & \Xi_{R,r} \\
499     \end{array} \right)\left( \begin{array}{l}
500 gezelter 3337 \mathbf{v}_{R}^b (t) \\
501     \mathbf{\omega} (t) \\
502 tim 2746 \end{array} \right),
503     \end{equation}
504 gezelter 3340 while the random force, $\mathbf{F}_{r}$, is a Gaussian stochastic
505     variable with zero mean and variance
506 tim 2746 \begin{equation}
507 gezelter 3340 \left\langle {\mathbf{F}_{r}(t) (\mathbf{F}_{r}(t'))^T } \right\rangle =
508 xsun 3338 \left\langle {\mathbf{F}_{r}^b (t) (\mathbf{F}_{r}^b (t'))^T } \right\rangle =
509 gezelter 3333 2 k_B T \Xi_R \delta(t - t'). \label{randomForce}
510 tim 2746 \end{equation}
511 gezelter 3340 $\Xi_R$ is the $6\times6$ resistance tensor at the center of
512     resistance. Once this tensor is known for a given rigid body,
513     obtaining a stochastic vector that has the properties in
514     Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a
515     one-time Cholesky decomposition to obtain the square root matrix of
516     the resistance tensor $\Xi_R = \mathbf{S} \mathbf{S}^{T}$, where
517     $\mathbf{S}$ is a lower triangular matrix.\cite{SchlickBook} A vector
518     with the statistics required for the random force can then be obtained
519     by multiplying $\mathbf{S}$ onto a 6-vector $Z$ which has elements
520     chosen from a Gaussian distribution, such that:
521     \begin{equation}
522     \langle Z_i \rangle = 0, \hspace{1in} \langle Z_i \cdot Z_j \rangle = \frac{2 k_B
523     T}{\delta t} \delta_{ij}.
524     \end{equation}
525     The random force, $F_{r}^{b} = \mathbf{S} Z$, can be shown to have the
526     correct ohmic
527 gezelter 3333
528 gezelter 3340
529     Each
530     time a random force vector is needed, a gaussian random vector is
531     generated and then the square root matrix is multiplied onto this
532     vector.
533    
534 gezelter 3337 The equation of motion for $\mathbf{v}$ can be written as
535 tim 2746 \begin{equation}
536 xsun 3339 m \dot{\mathbf{v}} (t) = \mathbf{f}_{s}^l (t) + \mathbf{f}_{f}^l (t) +
537 xsun 3338 \mathbf{f}_{r}^l (t)
538 tim 2746 \end{equation}
539     Since the frictional force is applied at the center of resistance
540     which generally does not coincide with the center of mass, an extra
541     torque is exerted at the center of mass. Thus, the net body-fixed
542 gezelter 3337 frictional torque at the center of mass, $\tau_{f}^b (t)$, is
543 tim 2746 given by
544     \begin{equation}
545 xsun 3339 \tau_{f}^b \leftarrow \tau_{f}^b + \mathbf{r}_{MR} \times \mathbf{f}_{f}^b
546 tim 2746 \end{equation}
547     where $r_{MR}$ is the vector from the center of mass to the center
548 tim 2999 of the resistance. Instead of integrating the angular velocity in
549     lab-fixed frame, we consider the equation of angular momentum in
550     body-fixed frame
551 tim 2746 \begin{equation}
552 xsun 3339 \dot j(t) = \tau_{s}^b (t) + \tau_{f}^b (t) + \tau_{r}^b(t)
553 tim 2746 \end{equation}
554 gezelter 3333 Embedding the friction terms into force and torque, one can integrate
555     the Langevin equations of motion for rigid body of arbitrary shape in
556     a velocity-Verlet style 2-part algorithm, where $h= \delta t$:
557 tim 2746
558 tim 2999 {\tt moveA:}
559 tim 2746 \begin{align*}
560 tim 2999 {\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t)
561     + \frac{h}{2} \left( {\bf f}(t) / m \right), \\
562     %
563     {\bf r}(t + h) &\leftarrow {\bf r}(t)
564     + h {\bf v}\left(t + h / 2 \right), \\
565     %
566     {\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t)
567     + \frac{h}{2} {\bf \tau}^b(t), \\
568     %
569     \mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
570     (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
571 tim 2746 \end{align*}
572     In this context, the $\mathrm{rotate}$ function is the reversible
573 tim 2999 product of the three body-fixed rotations,
574 tim 2746 \begin{equation}
575     \mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
576     \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
577     / 2) \cdot \mathsf{G}_x(a_x /2),
578     \end{equation}
579     where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
580 tim 2999 rotates both the rotation matrix ($\mathsf{Q}$) and the body-fixed
581     angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
582     axis $\alpha$,
583 tim 2746 \begin{equation}
584     \mathsf{G}_\alpha( \theta ) = \left\{
585     \begin{array}{lcl}
586 tim 2999 \mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
587 tim 2746 {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
588     j}(0).
589     \end{array}
590     \right.
591     \end{equation}
592     $\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
593     rotation matrix. For example, in the small-angle limit, the
594     rotation matrix around the body-fixed x-axis can be approximated as
595     \begin{equation}
596     \mathsf{R}_x(\theta) \approx \left(
597     \begin{array}{ccc}
598     1 & 0 & 0 \\
599     0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+
600     \theta^2 / 4} \\
601     0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
602     \theta^2 / 4}
603     \end{array}
604     \right).
605     \end{equation}
606 tim 2999 All other rotations follow in a straightforward manner. After the
607     first part of the propagation, the forces and body-fixed torques are
608     calculated at the new positions and orientations
609 tim 2746
610 tim 2999 {\tt doForces:}
611     \begin{align*}
612     {\bf f}(t + h) &\leftarrow
613     - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
614     %
615     {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
616     \times \frac{\partial V}{\partial {\bf u}}, \\
617     %
618     {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
619     \cdot {\bf \tau}^s(t + h).
620     \end{align*}
621 tim 2746 Once the forces and torques have been obtained at the new time step,
622     the velocities can be advanced to the same time value.
623    
624 tim 2999 {\tt moveB:}
625 tim 2746 \begin{align*}
626 tim 2999 {\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2
627     \right)
628     + \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\
629     %
630     {\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2
631     \right)
632     + \frac{h}{2} {\bf \tau}^b(t + h) .
633 tim 2746 \end{align*}
634    
635 gezelter 3310 \section{Validating the Method\label{sec:validating}}
636 gezelter 3302 In order to validate our Langevin integrator for arbitrarily-shaped
637 gezelter 3305 rigid bodies, we implemented the algorithm in {\sc
638     oopse}\cite{Meineke2005} and compared the results of this algorithm
639     with the known
640 gezelter 3302 hydrodynamic limiting behavior for a few model systems, and to
641     microcanonical molecular dynamics simulations for some more
642     complicated bodies. The model systems and their analytical behavior
643     (if known) are summarized below. Parameters for the primary particles
644     comprising our model systems are given in table \ref{tab:parameters},
645     and a sketch of the arrangement of these primary particles into the
646 gezelter 3305 model rigid bodies is shown in figure \ref{fig:models}. In table
647     \ref{tab:parameters}, $d$ and $l$ are the physical dimensions of
648     ellipsoidal (Gay-Berne) particles. For spherical particles, the value
649     of the Lennard-Jones $\sigma$ parameter is the particle diameter
650     ($d$). Gay-Berne ellipsoids have an energy scaling parameter,
651     $\epsilon^s$, which describes the well depth for two identical
652     ellipsoids in a {\it side-by-side} configuration. Additionally, a
653     well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$,
654     describes the ratio between the well depths in the {\it end-to-end}
655     and side-by-side configurations. For spheres, $\epsilon^r \equiv 1$.
656     Moments of inertia are also required to describe the motion of primary
657     particles with orientational degrees of freedom.
658 gezelter 3299
659 gezelter 3302 \begin{table*}
660     \begin{minipage}{\linewidth}
661     \begin{center}
662     \caption{Parameters for the primary particles in use by the rigid body
663     models in figure \ref{fig:models}.}
664     \begin{tabular}{lrcccccccc}
665     \hline
666     & & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\
667     & & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ (kcal/mol) & $\epsilon^r$ &
668     $m$ (amu) & $I_{xx}$ & $I_{yy}$ & $I_{zz}$ \\ \hline
669 gezelter 3308 Sphere & & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\
670 gezelter 3302 Ellipsoid & & 4.6 & 13.8 & 0.8 & 0.2 & 200 & 2105 & 2105 & 421 \\
671 gezelter 3308 Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\
672 gezelter 3302 Banana &(3 identical ellipsoids)& 4.2 & 11.2 & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\
673     Lipid: & Spherical Head & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\
674     & Ellipsoidal Tail & 4.6 & 13.8 & 0.8 & 0.2 & 760 & 45000 & 45000 & 9000 \\
675     Solvent & & 4.7 & $= d$ & 0.8 & 1 & 72.06 & & & \\
676     \hline
677     \end{tabular}
678     \label{tab:parameters}
679     \end{center}
680     \end{minipage}
681     \end{table*}
682    
683 gezelter 3305 \begin{figure}
684     \centering
685     \includegraphics[width=3in]{sketch}
686     \caption[Sketch of the model systems]{A sketch of the model systems
687     used in evaluating the behavior of the rigid body Langevin
688     integrator.} \label{fig:models}
689     \end{figure}
690    
691 gezelter 3302 \subsection{Simulation Methodology}
692     We performed reference microcanonical simulations with explicit
693     solvents for each of the different model system. In each case there
694     was one solute model and 1929 solvent molecules present in the
695     simulation box. All simulations were equilibrated using a
696     constant-pressure and temperature integrator with target values of 300
697     K for the temperature and 1 atm for pressure. Following this stage,
698     further equilibration and sampling was done in a microcanonical
699 gezelter 3305 ensemble. Since the model bodies are typically quite massive, we were
700 gezelter 3310 able to use a time step of 25 fs.
701    
702     The model systems studied used both Lennard-Jones spheres as well as
703     uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne
704     potential was a single site model for the interactions of rigid
705     ellipsoidal molecules.\cite{Gay81} It can be thought of as a
706     modification of the Gaussian overlap model originally described by
707     Berne and Pechukas.\cite{Berne72} The potential is constructed in the
708     familiar form of the Lennard-Jones function using
709     orientation-dependent $\sigma$ and $\epsilon$ parameters,
710     \begin{equation*}
711     V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
712     r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
713     {\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u
714     }_i},
715     {\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
716     -\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
717     {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
718     \label{eq:gb}
719     \end{equation*}
720    
721     The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
722     \hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf
723     \hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters
724     are dependent on the relative orientations of the two ellipsoids (${\bf
725     \hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the
726     inter-ellipsoid separation (${\bf \hat{r}}_{ij}$). The shape and
727     attractiveness of each ellipsoid is governed by a relatively small set
728     of parameters: $l$ and $d$ describe the length and width of each
729     uniaxial ellipsoid, while $\epsilon^s$, which describes the well depth
730     for two identical ellipsoids in a {\it side-by-side} configuration.
731     Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e /
732     \epsilon^s$, describes the ratio between the well depths in the {\it
733     end-to-end} and side-by-side configurations. Details of the potential
734     are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an
735     excellent overview of the computational methods that can be used to
736     efficiently compute forces and torques for this potential can be found
737     in Ref. \citen{Golubkov06}
738    
739     For the interaction between nonequivalent uniaxial ellipsoids (or
740     between spheres and ellipsoids), the spheres are treated as ellipsoids
741     with an aspect ratio of 1 ($d = l$) and with an well depth ratio
742     ($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$). The form of the
743     Gay-Berne potential we are using was generalized by Cleaver {\it et
744     al.} and is appropriate for dissimilar uniaxial
745     ellipsoids.\cite{Cleaver96}
746    
747     A switching function was applied to all potentials to smoothly turn
748     off the interactions between a range of $22$ and $25$ \AA. The
749     switching function was the standard (cubic) function,
750 gezelter 3302 \begin{equation}
751     s(r) =
752     \begin{cases}
753     1 & \text{if $r \le r_{\text{sw}}$},\\
754     \frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2}
755     {(r_{\text{cut}} - r_{\text{sw}})^3}
756     & \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\
757     0 & \text{if $r > r_{\text{cut}}$.}
758     \end{cases}
759     \label{eq:switchingFunc}
760     \end{equation}
761 gezelter 3310
762 gezelter 3302 To measure shear viscosities from our microcanonical simulations, we
763     used the Einstein form of the pressure correlation function,\cite{hess:209}
764     \begin{equation}
765 gezelter 3310 \eta = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left(
766     \int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}.
767 gezelter 3302 \label{eq:shear}
768     \end{equation}
769     A similar form exists for the bulk viscosity
770     \begin{equation}
771 gezelter 3310 \kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left(
772 gezelter 3302 \int_{t_0}^{t_0 + t}
773 gezelter 3310 \left(P\left(t'\right)-\left\langle P \right\rangle \right)dt'
774     \right)^2 \right\rangle_{t_0}.
775 gezelter 3302 \end{equation}
776     Alternatively, the shear viscosity can also be calculated using a
777     Green-Kubo formula with the off-diagonal pressure tensor correlation function,
778     \begin{equation}
779 gezelter 3310 \eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0
780     + t) \right\rangle_{t_0} dt,
781 gezelter 3302 \end{equation}
782     although this method converges extremely slowly and is not practical
783     for obtaining viscosities from molecular dynamics simulations.
784    
785     The Langevin dynamics for the different model systems were performed
786     at the same temperature as the average temperature of the
787     microcanonical simulations and with a solvent viscosity taken from
788 gezelter 3305 Eq. (\ref{eq:shear}) applied to these simulations. We used 1024
789     independent solute simulations to obtain statistics on our Langevin
790     integrator.
791 gezelter 3302
792     \subsection{Analysis}
793    
794     The quantities of interest when comparing the Langevin integrator to
795     analytic hydrodynamic equations and to molecular dynamics simulations
796     are typically translational diffusion constants and orientational
797     relaxation times. Translational diffusion constants for point
798     particles are computed easily from the long-time slope of the
799     mean-square displacement,
800     \begin{equation}
801 gezelter 3310 D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \left\langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \right\rangle,
802 gezelter 3302 \end{equation}
803     of the solute molecules. For models in which the translational
804 gezelter 3305 diffusion tensor (${\bf D}_{tt}$) has non-degenerate eigenvalues
805     (i.e. any non-spherically-symmetric rigid body), it is possible to
806     compute the diffusive behavior for motion parallel to each body-fixed
807     axis by projecting the displacement of the particle onto the
808     body-fixed reference frame at $t=0$. With an isotropic solvent, as we
809     have used in this study, there are differences between the three
810 gezelter 3302 diffusion constants, but these must converge to the same value at
811     longer times. Translational diffusion constants for the different
812 gezelter 3305 shaped models are shown in table \ref{tab:translation}.
813 gezelter 3302
814 gezelter 3305 In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational
815 gezelter 3302 diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object
816     {\it around} a particular body-fixed axis and {\it not} the diffusion
817     of a vector pointing along the axis. However, these eigenvalues can
818     be combined to find 5 characteristic rotational relaxation
819 gezelter 3305 times,\cite{PhysRev.119.53,Berne90}
820 gezelter 3302 \begin{eqnarray}
821 gezelter 3305 1 / \tau_1 & = & 6 D_r + 2 \Delta \\
822     1 / \tau_2 & = & 6 D_r - 2 \Delta \\
823     1 / \tau_3 & = & 3 (D_r + D_1) \\
824     1 / \tau_4 & = & 3 (D_r + D_2) \\
825     1 / \tau_5 & = & 3 (D_r + D_3)
826 gezelter 3302 \end{eqnarray}
827     where
828     \begin{equation}
829     D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right)
830     \end{equation}
831     and
832     \begin{equation}
833 gezelter 3305 \Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2}
834 gezelter 3302 \end{equation}
835 gezelter 3305 Each of these characteristic times can be used to predict the decay of
836     part of the rotational correlation function when $\ell = 2$,
837 gezelter 3302 \begin{equation}
838 gezelter 3305 C_2(t) = \frac{a^2}{N^2} e^{-t/\tau_1} + \frac{b^2}{N^2} e^{-t/\tau_2}.
839 gezelter 3302 \end{equation}
840 gezelter 3305 This is the same as the $F^2_{0,0}(t)$ correlation function that
841     appears in Ref. \citen{Berne90}. The amplitudes of the two decay
842     terms are expressed in terms of three dimensionless functions of the
843     eigenvalues: $a = \sqrt{3} (D_1 - D_2)$, $b = (2D_3 - D_1 - D_2 +
844     2\Delta)$, and $N = 2 \sqrt{\Delta b}$. Similar expressions can be
845     obtained for other angular momentum correlation
846     functions.\cite{PhysRev.119.53,Berne90} In all of the model systems we
847     studied, only one of the amplitudes of the two decay terms was
848     non-zero, so it was possible to derive a single relaxation time for
849     each of the hydrodynamic tensors. In many cases, these characteristic
850     times are averaged and reported in the literature as a single relaxation
851     time,\cite{Garcia-de-la-Torre:1997qy}
852 gezelter 3302 \begin{equation}
853 gezelter 3305 1 / \tau_0 = \frac{1}{5} \sum_{i=1}^5 \tau_{i}^{-1},
854     \end{equation}
855     although for the cases reported here, this averaging is not necessary
856     and only one of the five relaxation times is relevant.
857    
858     To test the Langevin integrator's behavior for rotational relaxation,
859     we have compared the analytical orientational relaxation times (if
860     they are known) with the general result from the diffusion tensor and
861     with the results from both the explicitly solvated molecular dynamics
862     and Langevin simulations. Relaxation times from simulations (both
863     microcanonical and Langevin), were computed using Legendre polynomial
864     correlation functions for a unit vector (${\bf u}$) fixed along one or
865     more of the body-fixed axes of the model.
866     \begin{equation}
867 gezelter 3310 C_{\ell}(t) = \left\langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf
868     u}_{i}(0) \right) \right\rangle
869 gezelter 3302 \end{equation}
870     For simulations in the high-friction limit, orientational correlation
871     times can then be obtained from exponential fits of this function, or by
872     integrating,
873     \begin{equation}
874 gezelter 3305 \tau = \ell (\ell + 1) \int_0^{\infty} C_{\ell}(t) dt.
875 gezelter 3302 \end{equation}
876 gezelter 3305 In lower-friction solvents, the Legendre correlation functions often
877     exhibit non-exponential decay, and may not be characterized by a
878     single decay constant.
879 gezelter 3302
880     In table \ref{tab:rotation} we show the characteristic rotational
881     relaxation times (based on the diffusion tensor) for each of the model
882     systems compared with the values obtained via microcanonical and Langevin
883     simulations.
884    
885 gezelter 3305 \subsection{Spherical particles}
886 gezelter 3299 Our model system for spherical particles was a Lennard-Jones sphere of
887     diameter ($\sigma$) 6.5 \AA\ in a sea of smaller spheres ($\sigma$ =
888     4.7 \AA). The well depth ($\epsilon$) for both particles was set to
889 gezelter 3302 an arbitrary value of 0.8 kcal/mol.
890 gezelter 3299
891     The Stokes-Einstein behavior of large spherical particles in
892     hydrodynamic flows is well known, giving translational friction
893     coefficients of $6 \pi \eta R$ (stick boundary conditions) and
894 gezelter 3302 rotational friction coefficients of $8 \pi \eta R^3$. Recently,
895     Schmidt and Skinner have computed the behavior of spherical tag
896     particles in molecular dynamics simulations, and have shown that {\it
897     slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more
898 gezelter 3299 appropriate for molecule-sized spheres embedded in a sea of spherical
899 gezelter 3310 solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx}
900 gezelter 3299
901     Our simulation results show similar behavior to the behavior observed
902 gezelter 3302 by Schmidt and Skinner. The diffusion constant obtained from our
903 gezelter 3299 microcanonical molecular dynamics simulations lies between the slip
904     and stick boundary condition results obtained via Stokes-Einstein
905     behavior. Since the Langevin integrator assumes Stokes-Einstein stick
906     boundary conditions in calculating the drag and random forces for
907     spherical particles, our Langevin routine obtains nearly quantitative
908     agreement with the hydrodynamic results for spherical particles. One
909     avenue for improvement of the method would be to compute elements of
910     $\Xi_{tt}$ assuming behavior intermediate between the two boundary
911 gezelter 3302 conditions.
912 gezelter 3299
913 gezelter 3310 In the explicit solvent simulations, both our solute and solvent
914     particles were structureless, exerting no torques upon each other.
915     Therefore, there are not rotational correlation times available for
916     this model system.
917 gezelter 3299
918 gezelter 3310 \subsection{Ellipsoids}
919     For uniaxial ellipsoids ($a > b = c$), Perrin's formulae for both
920 gezelter 3299 translational and rotational diffusion of each of the body-fixed axes
921     can be combined to give a single translational diffusion
922 gezelter 3302 constant,\cite{Berne90}
923 gezelter 3299 \begin{equation}
924     D = \frac{k_B T}{6 \pi \eta a} G(\rho),
925     \label{Dperrin}
926     \end{equation}
927     as well as a single rotational diffusion coefficient,
928     \begin{equation}
929     \Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2)
930     G(\rho) - 1}{1 - \rho^4} \right\}.
931     \label{ThetaPerrin}
932     \end{equation}
933     In these expressions, $G(\rho)$ is a function of the axial ratio
934     ($\rho = b / a$), which for prolate ellipsoids, is
935     \begin{equation}
936     G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 -
937     \rho^2)^{1/2}}{\rho} \right\}
938     \label{GPerrin}
939     \end{equation}
940     Again, there is some uncertainty about the correct boundary conditions
941     to use for molecular-scale ellipsoids in a sea of similarly-sized
942     solvent particles. Ravichandran and Bagchi found that {\it slip}
943 gezelter 3302 boundary conditions most closely resembled the simulation
944     results,\cite{Ravichandran:1999fk} in agreement with earlier work of
945     Tang and Evans.\cite{TANG:1993lr}
946 gezelter 3299
947 gezelter 3305 Even though there are analytic resistance tensors for ellipsoids, we
948     constructed a rough-shell model using 2135 beads (each with a diameter
949 gezelter 3310 of 0.25 \AA) to approximate the shape of the model ellipsoid. We
950 gezelter 3305 compared the Langevin dynamics from both the simple ellipsoidal
951     resistance tensor and the rough shell approximation with
952     microcanonical simulations and the predictions of Perrin. As in the
953     case of our spherical model system, the Langevin integrator reproduces
954     almost exactly the behavior of the Perrin formulae (which is
955     unsurprising given that the Perrin formulae were used to derive the
956 gezelter 3299 drag and random forces applied to the ellipsoid). We obtain
957     translational diffusion constants and rotational correlation times
958     that are within a few percent of the analytic values for both the
959     exact treatment of the diffusion tensor as well as the rough-shell
960     model for the ellipsoid.
961    
962 gezelter 3308 The translational diffusion constants from the microcanonical simulations
963     agree well with the predictions of the Perrin model, although the rotational
964     correlation times are a factor of 2 shorter than expected from hydrodynamic
965     theory. One explanation for the slower rotation
966     of explicitly-solvated ellipsoids is the possibility that solute-solvent
967     collisions happen at both ends of the solute whenever the principal
968     axis of the ellipsoid is turning. In the upper portion of figure
969     \ref{fig:explanation} we sketch a physical picture of this explanation.
970     Since our Langevin integrator is providing nearly quantitative agreement with
971     the Perrin model, it also predicts orientational diffusion for ellipsoids that
972     exceed explicitly solvated correlation times by a factor of two.
973 gezelter 3299
974 gezelter 3310 \subsection{Rigid dumbbells}
975 gezelter 3302 Perhaps the only {\it composite} rigid body for which analytic
976     expressions for the hydrodynamic tensor are available is the
977     two-sphere dumbbell model. This model consists of two non-overlapping
978     spheres held by a rigid bond connecting their centers. There are
979     competing expressions for the 6x6 resistance tensor for this
980     model. Equation (\ref{introEquation:oseenTensor}) above gives the
981     original Oseen tensor, while the second order expression introduced by
982     Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la
983     Torre and Bloomfield,\cite{Torre1977} is given above as
984 gezelter 3299 Eq. (\ref{introEquation:RPTensorNonOverlapped}). In our case, we use
985     a model dumbbell in which the two spheres are identical Lennard-Jones
986     particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at
987 gezelter 3302 a distance of 6.532 \AA.
988 gezelter 3299
989     The theoretical values for the translational diffusion constant of the
990     dumbbell are calculated from the work of Stimson and Jeffery, who
991     studied the motion of this system in a flow parallel to the
992 gezelter 3302 inter-sphere axis,\cite{Stimson:1926qy} and Davis, who studied the
993     motion in a flow {\it perpendicular} to the inter-sphere
994     axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it
995     orientational} correlation times for this model system (other than
996 gezelter 3305 those derived from the 6 x 6 tensors mentioned above).
997 gezelter 3299
998 gezelter 3305 The bead model for this model system comprises the two large spheres
999     by themselves, while the rough shell approximation used 3368 separate
1000     beads (each with a diameter of 0.25 \AA) to approximate the shape of
1001     the rigid body. The hydrodynamics tensors computed from both the bead
1002     and rough shell models are remarkably similar. Computing the initial
1003     hydrodynamic tensor for a rough shell model can be quite expensive (in
1004     this case it requires inverting a 10104 x 10104 matrix), while the
1005     bead model is typically easy to compute (in this case requiring
1006 gezelter 3308 inversion of a 6 x 6 matrix).
1007 gezelter 3305
1008 gezelter 3308 \begin{figure}
1009     \centering
1010 gezelter 3310 \includegraphics[width=2in]{RoughShell}
1011 gezelter 3308 \caption[Model rigid bodies and their rough shell approximations]{The
1012     model rigid bodies (left column) used to test this algorithm and their
1013     rough-shell approximations (right-column) that were used to compute
1014     the hydrodynamic tensors. The top two models (ellipsoid and dumbbell)
1015     have analytic solutions and were used to test the rough shell
1016     approximation. The lower two models (banana and lipid) were compared
1017     with explicitly-solvated molecular dynamics simulations. }
1018     \label{fig:roughShell}
1019     \end{figure}
1020    
1021    
1022 gezelter 3305 Once the hydrodynamic tensor has been computed, there is no additional
1023     penalty for carrying out a Langevin simulation with either of the two
1024     different hydrodynamics models. Our naive expectation is that since
1025     the rigid body's surface is roughened under the various shell models,
1026     the diffusion constants will be even farther from the ``slip''
1027     boundary conditions than observed for the bead model (which uses a
1028     Stokes-Einstein model to arrive at the hydrodynamic tensor). For the
1029     dumbbell, this prediction is correct although all of the Langevin
1030     diffusion constants are within 6\% of the diffusion constant predicted
1031     from the fully solvated system.
1032    
1033 gezelter 3308 For rotational motion, Langevin integration (and the hydrodynamic tensor)
1034     yields rotational correlation times that are substantially shorter than those
1035     obtained from explicitly-solvated simulations. It is likely that this is due
1036     to the large size of the explicit solvent spheres, a feature that prevents
1037     the solvent from coming in contact with a substantial fraction of the surface
1038     area of the dumbbell. Therefore, the explicit solvent only provides drag
1039     over a substantially reduced surface area of this model, while the
1040     hydrodynamic theories utilize the entire surface area for estimating
1041     rotational diffusion. A sketch of the free volume available in the explicit
1042     solvent simulations is shown in figure \ref{fig:explanation}.
1043 gezelter 3305
1044 gezelter 3310
1045     \begin{figure}
1046     \centering
1047     \includegraphics[width=6in]{explanation}
1048     \caption[Explanations of the differences between orientational
1049     correlation times for explicitly-solvated models and hydrodynamics
1050     predictions]{Explanations of the differences between orientational
1051     correlation times for explicitly-solvated models and hydrodynamic
1052     predictions. For the ellipsoids (upper figures), rotation of the
1053     principal axis can involve correlated collisions at both sides of the
1054     solute. In the rigid dumbbell model (lower figures), the large size
1055     of the explicit solvent spheres prevents them from coming in contact
1056     with a substantial fraction of the surface area of the dumbbell.
1057     Therefore, the explicit solvent only provides drag over a
1058     substantially reduced surface area of this model, where the
1059     hydrodynamic theories utilize the entire surface area for estimating
1060     rotational diffusion.
1061     } \label{fig:explanation}
1062     \end{figure}
1063    
1064    
1065    
1066     \subsection{Composite banana-shaped molecules}
1067     Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids have
1068     been used by Orlandi {\it et al.} to observe mesophases in
1069     coarse-grained models for bent-core liquid crystalline
1070     molecules.\cite{Orlandi:2006fk} We have used the same overlapping
1071 gezelter 3299 ellipsoids as a way to test the behavior of our algorithm for a
1072     structure of some interest to the materials science community,
1073     although since we are interested in capturing only the hydrodynamic
1074 gezelter 3310 behavior of this model, we have left out the dipolar interactions of
1075     the original Orlandi model.
1076 gezelter 3308
1077     A reference system composed of a single banana rigid body embedded in a
1078     sea of 1929 solvent particles was created and run under standard
1079     (microcanonical) molecular dynamics. The resulting viscosity of this
1080     mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})).
1081     To calculate the hydrodynamic properties of the banana rigid body model,
1082 gezelter 3310 we created a rough shell (see Fig.~\ref{fig:roughShell}), in which
1083 gezelter 3308 the banana is represented as a ``shell'' made of 3321 identical beads
1084 gezelter 3310 (0.25 \AA\ in diameter) distributed on the surface. Applying the
1085 gezelter 3308 procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
1086 gezelter 3310 identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 \AA), as
1087     well as the resistance tensor,
1088     \begin{equation*}
1089     \Xi =
1090 gezelter 3308 \left( {\begin{array}{*{20}c}
1091     0.9261 & 0 & 0&0&0.08585&0.2057\\
1092     0& 0.9270&-0.007063& 0.08585&0&0\\
1093     0&-0.007063&0.7494&0.2057&0&0\\
1094 gezelter 3310 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),
1095     \end{equation*}
1096     where the units for translational, translation-rotation coupling and
1097     rotational tensors are (kcal fs / mol \AA$^2$), (kcal fs / mol \AA\ rad),
1098     and (kcal fs / mol rad$^2$), respectively.
1099 gezelter 3299
1100 gezelter 3308 The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor)
1101     are essentially quantitative for translational diffusion of this model.
1102     Orientational correlation times under the Langevin rigid-body integrator
1103     are within 11\% of the values obtained from explicit solvent, but these
1104     models also exhibit some solvent inaccessible surface area in the
1105     explicitly-solvated case.
1106    
1107 gezelter 3310 \subsection{Composite sphero-ellipsoids}
1108 gezelter 3299 Spherical heads perched on the ends of Gay-Berne ellipsoids have been
1109 xsun 3312 used recently as models for lipid
1110     molecules.\cite{SunGezelter08,Ayton01}
1111 gezelter 3310 MORE DETAILS
1112 xsun 3298
1113 xsun 3312 A reference system composed of a single lipid rigid body embedded in a
1114     sea of 1929 solvent particles was created and run under standard
1115     (microcanonical) molecular dynamics. The resulting viscosity of this
1116     mixture was 0.349 centipoise (as estimated using
1117     Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties of
1118     the lipid rigid body model, we created a rough shell (see
1119     Fig.~\ref{fig:roughShell}), in which the lipid is represented as a
1120     ``shell'' made of 3550 identical beads (0.25 \AA\ in diameter)
1121     distributed on the surface. Applying the procedure described in
1122     Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
1123     identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46
1124     \AA).
1125 gezelter 3310
1126 gezelter 3315
1127 gezelter 3310 \subsection{Summary}
1128 xsun 3298 According to our simulations, the langevin dynamics is a reliable
1129     theory to apply to replace the explicit solvents, especially for the
1130     translation properties. For large molecules, the rotation properties
1131     are also mimiced reasonablly well.
1132    
1133 gezelter 3315 \begin{figure}
1134     \centering
1135     \includegraphics[width=\linewidth]{graph}
1136     \caption[Mean squared displacements and orientational
1137     correlation functions for each of the model rigid bodies.]{The
1138     mean-squared displacements ($\langle r^2(t) \rangle$) and
1139     orientational correlation functions ($C_2(t)$) for each of the model
1140     rigid bodies studied. The circles are the results for microcanonical
1141     simulations with explicit solvent molecules, while the other data sets
1142     are results for Langevin dynamics using the different hydrodynamic
1143     tensor approximations. The Perrin model for the ellipsoids is
1144     considered the ``exact'' hydrodynamic behavior (this can also be said
1145     for the translational motion of the dumbbell operating under the bead
1146     model). In most cases, the various hydrodynamics models reproduce
1147     each other quantitatively.}
1148     \label{fig:results}
1149     \end{figure}
1150    
1151 xsun 3298 \begin{table*}
1152     \begin{minipage}{\linewidth}
1153     \begin{center}
1154 gezelter 3305 \caption{Translational diffusion constants (D) for the model systems
1155     calculated using microcanonical simulations (with explicit solvent),
1156     theoretical predictions, and Langevin simulations (with implicit solvent).
1157     Analytical solutions for the exactly-solved hydrodynamics models are
1158     from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936}
1159     (ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq}
1160     (dumbbell). The other model systems have no known analytic solution.
1161     All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (=
1162     $10^{-4}$ \AA$^2$ / fs). }
1163     \begin{tabular}{lccccccc}
1164 xsun 3298 \hline
1165 gezelter 3305 & \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\
1166     \cline{2-3} \cline{5-7}
1167     model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\
1168 xsun 3298 \hline
1169 xsun 3312 sphere & 0.279 & 3.06 & & 2.42 & exact & 2.42 & 2.33 \\
1170 gezelter 3305 ellipsoid & 0.255 & 2.44 & & 2.34 & exact & 2.34 & 2.37 \\
1171     & 0.255 & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\
1172 xsun 3312 dumbbell & 0.308 & 2.06 & & 1.64 & bead model & 1.65 & 1.62 \\
1173     & 0.308 & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\
1174 gezelter 3305 banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\
1175     lipid & 0.349 & 0.96 & & & rough shell & 1.33 & 1.32 \\
1176 xsun 3298 \end{tabular}
1177     \label{tab:translation}
1178     \end{center}
1179     \end{minipage}
1180     \end{table*}
1181    
1182     \begin{table*}
1183     \begin{minipage}{\linewidth}
1184     \begin{center}
1185 gezelter 3305 \caption{Orientational relaxation times ($\tau$) for the model systems using
1186     microcanonical simulation (with explicit solvent), theoretical
1187     predictions, and Langevin simulations (with implicit solvent). All
1188     relaxation times are for the rotational correlation function with
1189     $\ell = 2$ and are reported in units of ps. The ellipsoidal model has
1190     an exact solution for the orientational correlation time due to
1191     Perrin, but the other model systems have no known analytic solution.}
1192     \begin{tabular}{lccccccc}
1193 xsun 3298 \hline
1194 gezelter 3305 & \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\
1195     \cline{2-3} \cline{5-7}
1196     model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic & simulation \\
1197 xsun 3298 \hline
1198 xsun 3312 sphere & 0.279 & & & 9.69 & exact & 9.69 & 9.64 \\
1199 gezelter 3305 ellipsoid & 0.255 & 46.7 & & 22.0 & exact & 22.0 & 22.2 \\
1200     & 0.255 & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\
1201 xsun 3312 dumbbell & 0.308 & 14.1 & & & bead model & 50.0 & 50.1 \\
1202     & 0.308 & 14.1 & & & rough shell & 41.5 & 41.3 \\
1203 gezelter 3305 banana & 0.298 & 63.8 & & & rough shell & 70.9 & 70.9 \\
1204     lipid & 0.349 & 78.0 & & & rough shell & 76.9 & 77.9 \\
1205     \hline
1206 xsun 3298 \end{tabular}
1207     \label{tab:rotation}
1208     \end{center}
1209     \end{minipage}
1210     \end{table*}
1211    
1212 gezelter 3310 \section{Application: A rigid-body lipid bilayer}
1213    
1214     The Langevin dynamics integrator was applied to study the formation of
1215     corrugated structures emerging from simulations of the coarse grained
1216     lipid molecular models presented above. The initial configuration is
1217 xsun 3298 taken from our molecular dynamics studies on lipid bilayers with
1218 gezelter 3310 lennard-Jones sphere solvents. The solvent molecules were excluded
1219     from the system, and the experimental value for the viscosity of water
1220     at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects
1221     of the solvent. The absence of explicit solvent molecules and the
1222     stability of the integrator allowed us to take timesteps of 50 fs. A
1223     total simulation run time of 100 ns was sampled.
1224     Fig. \ref{fig:bilayer} shows the configuration of the system after 100
1225     ns, and the ripple structure remains stable during the entire
1226     trajectory. Compared with using explicit bead-model solvent
1227     molecules, the efficiency of the simulation has increased by an order
1228 xsun 3298 of magnitude.
1229    
1230 gezelter 3310 \begin{figure}
1231     \centering
1232     \includegraphics[width=\linewidth]{bilayer}
1233     \caption[Snapshot of a bilayer of rigid-body models for lipids]{A
1234     snapshot of a bilayer composed of rigid-body models for lipid
1235     molecules evolving using the Langevin integrator described in this
1236     work.} \label{fig:bilayer}
1237     \end{figure}
1238    
1239 tim 2746 \section{Conclusions}
1240    
1241 tim 2999 We have presented a new Langevin algorithm by incorporating the
1242     hydrodynamics properties of arbitrary shaped molecules into an
1243 gezelter 3308 advanced symplectic integration scheme. Further studies in systems
1244     involving banana shaped molecules illustrated that the dynamic
1245     properties could be preserved by using this new algorithm as an
1246     implicit solvent model.
1247 tim 2999
1248    
1249 tim 2746 \section{Acknowledgments}
1250     Support for this project was provided by the National Science
1251     Foundation under grant CHE-0134881. T.L. also acknowledges the
1252     financial support from center of applied mathematics at University
1253     of Notre Dame.
1254     \newpage
1255    
1256 gezelter 3305 \bibliographystyle{jcp}
1257 tim 2746 \bibliography{langevin}
1258    
1259     \end{document}