1 |
tim |
2746 |
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
2 |
|
|
%\documentclass[aps,prb,preprint]{revtex4} |
3 |
|
|
\documentclass[11pt]{article} |
4 |
|
|
\usepackage{endfloat} |
5 |
|
|
\usepackage{amsmath,bm} |
6 |
|
|
\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 |
tim |
2999 |
As alternative to Newtonian dynamics, Langevin dynamics, which |
52 |
|
|
mimics a simple heat bath with stochastic and dissipative forces, |
53 |
|
|
has been applied in a variety of studies. The stochastic treatment |
54 |
|
|
of the solvent enables us to carry out substantially longer time |
55 |
|
|
simulations. Implicit solvent Langevin dynamics simulations of |
56 |
|
|
met-enkephalin not only outperform explicit solvent simulations for |
57 |
|
|
computational efficiency, but also agrees very well with explicit |
58 |
|
|
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 |
tim |
2746 |
instance, the oscillation power spectrum of nanoparticles from |
76 |
|
|
Langevin dynamics simulation has the same peak frequencies for |
77 |
tim |
2999 |
different wave vectors, which recovers the property of magnetic |
78 |
|
|
excitations in small finite structures.\cite{Berkov2005a} |
79 |
tim |
2746 |
|
80 |
|
|
%review rigid body dynamics |
81 |
|
|
Rigid bodies are frequently involved in the modeling of different |
82 |
|
|
areas, from engineering, physics, to chemistry. For example, |
83 |
|
|
missiles and vehicle are usually modeled by rigid bodies. The |
84 |
|
|
movement of the objects in 3D gaming engine or other physics |
85 |
|
|
simulator is governed by the rigid body dynamics. In molecular |
86 |
|
|
simulation, rigid body is used to simplify the model in |
87 |
|
|
protein-protein docking study{\cite{Gray2003}}. |
88 |
|
|
|
89 |
|
|
It is very important to develop stable and efficient methods to |
90 |
tim |
2999 |
integrate the equations of motion for orientational degrees of |
91 |
|
|
freedom. Euler angles are the natural choice to describe the |
92 |
|
|
rotational degrees of freedom. However, due to $\frac {1}{sin |
93 |
|
|
\theta}$ singularities, the numerical integration of corresponding |
94 |
|
|
equations of these motion is very inefficient and inaccurate. |
95 |
|
|
Although an alternative integrator using multiple sets of Euler |
96 |
|
|
angles can overcome this difficulty\cite{Barojas1973}, the |
97 |
|
|
computational penalty and the loss of angular momentum conservation |
98 |
|
|
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} |
109 |
tim |
2746 |
|
110 |
tim |
2999 |
A break-through in geometric literature suggests that, in order to |
111 |
tim |
2746 |
develop a long-term integration scheme, one should preserve the |
112 |
tim |
2999 |
symplectic structure of the propagator. By introducing a conjugate |
113 |
|
|
momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's |
114 |
|
|
equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was |
115 |
|
|
proposed to evolve the Hamiltonian system in a constraint manifold |
116 |
|
|
by iteratively satisfying the orthogonality constraint $Q^T Q = 1$. |
117 |
|
|
An alternative method using the quaternion representation was |
118 |
|
|
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. |
122 |
tim |
2746 |
|
123 |
|
|
%review langevin/browninan dynamics for arbitrarily shaped rigid body |
124 |
|
|
Combining Langevin or Brownian dynamics with rigid body dynamics, |
125 |
tim |
2999 |
one can study slow processes in biomolecular systems. Modeling DNA |
126 |
|
|
as a chain of rigid beads, which are subject to harmonic potentials |
127 |
|
|
as well as excluded volume potentials, Mielke and his coworkers |
128 |
|
|
discovered rapid superhelical stress generations from the stochastic |
129 |
|
|
simulation of twin supercoiling DNA with response to induced |
130 |
|
|
torques.\cite{Mielke2004} Membrane fusion is another key biological |
131 |
|
|
process which controls a variety of physiological functions, such as |
132 |
|
|
release of neurotransmitters \textit{etc}. A typical fusion event |
133 |
|
|
happens on the time scale of a millisecond, which is impractical to |
134 |
|
|
study using atomistic models with newtonian mechanics. With the help |
135 |
|
|
of coarse-grained rigid body model and stochastic dynamics, the |
136 |
|
|
fusion pathways were explored by many |
137 |
|
|
researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the |
138 |
|
|
difficulty of numerical integration of anisotropic rotation, most of |
139 |
|
|
the rigid body models are simply modeled using spheres, cylinders, |
140 |
|
|
ellipsoids or other regular shapes in stochastic simulations. In an |
141 |
|
|
effort to account for the diffusion anisotropy of arbitrary |
142 |
tim |
2746 |
particles, Fernandes and de la Torre improved the original Brownian |
143 |
|
|
dynamics simulation algorithm\cite{Ermak1978,Allison1991} by |
144 |
|
|
incorporating a generalized $6\times6$ diffusion tensor and |
145 |
|
|
introducing a simple rotation evolution scheme consisting of three |
146 |
tim |
2999 |
consecutive rotations.\cite{Fernandes2002} Unfortunately, unexpected |
147 |
|
|
errors and biases are introduced into the system due to the |
148 |
|
|
arbitrary order of applying the noncommuting rotation |
149 |
|
|
operators.\cite{Beard2003} Based on the observation the momentum |
150 |
tim |
2746 |
relaxation time is much less than the time step, one may ignore the |
151 |
tim |
2999 |
inertia in Brownian dynamics. However, the assumption of zero |
152 |
tim |
2746 |
average acceleration is not always true for cooperative motion which |
153 |
|
|
is common in protein motion. An inertial Brownian dynamics (IBD) was |
154 |
|
|
proposed to address this issue by adding an inertial correction |
155 |
tim |
2999 |
term.\cite{Beard2000} As a complement to IBD which has a lower bound |
156 |
tim |
2746 |
in time step because of the inertial relaxation time, long-time-step |
157 |
|
|
inertial dynamics (LTID) can be used to investigate the inertial |
158 |
|
|
behavior of the polymer segments in low friction |
159 |
tim |
2999 |
regime.\cite{Beard2000} LTID can also deal with the rotational |
160 |
tim |
2746 |
dynamics for nonskew bodies without translation-rotation coupling by |
161 |
|
|
separating the translation and rotation motion and taking advantage |
162 |
|
|
of the analytical solution of hydrodynamics properties. However, |
163 |
tim |
2999 |
typical nonskew bodies like cylinders and ellipsoids are inadequate |
164 |
|
|
to represent most complex macromolecule assemblies. These intricate |
165 |
tim |
2746 |
molecules have been represented by a set of beads and their |
166 |
tim |
2999 |
hydrodynamic properties can be calculated using variants on the |
167 |
|
|
standard hydrodynamic interaction tensors. |
168 |
tim |
2746 |
|
169 |
|
|
The goal of the present work is to develop a Langevin dynamics |
170 |
tim |
2999 |
algorithm for arbitrary-shaped rigid particles by integrating the |
171 |
|
|
accurate estimation of friction tensor from hydrodynamics theory |
172 |
|
|
into the sophisticated rigid body dynamics algorithms. |
173 |
tim |
2746 |
|
174 |
tim |
2999 |
\section{Computational Methods{\label{methodSec}}} |
175 |
tim |
2746 |
|
176 |
tim |
2999 |
\subsection{\label{introSection:frictionTensor}Friction Tensor} |
177 |
|
|
Theoretically, the friction kernel can be determined using the |
178 |
|
|
velocity autocorrelation function. However, this approach becomes |
179 |
|
|
impractical when the system becomes more and more complicated. |
180 |
|
|
Instead, various approaches based on hydrodynamics have been |
181 |
|
|
developed to calculate the friction coefficients. In general, the |
182 |
|
|
friction tensor $\Xi$ is a $6\times 6$ matrix given by |
183 |
|
|
\[ |
184 |
|
|
\Xi = \left( {\begin{array}{*{20}c} |
185 |
|
|
{\Xi _{}^{tt} } & {\Xi _{}^{rt} } \\ |
186 |
|
|
{\Xi _{}^{tr} } & {\Xi _{}^{rr} } \\ |
187 |
|
|
\end{array}} \right). |
188 |
|
|
\] |
189 |
|
|
Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$ |
190 |
|
|
translational friction tensor and rotational resistance (friction) |
191 |
|
|
tensor respectively, while ${\Xi^{tr} }$ is translation-rotation |
192 |
|
|
coupling tensor and $ {\Xi^{rt} }$ is rotation-translation coupling |
193 |
|
|
tensor. When a particle moves in a fluid, it may experience friction |
194 |
|
|
force or torque along the opposite direction of the velocity or |
195 |
|
|
angular velocity, |
196 |
|
|
\[ |
197 |
tim |
2746 |
\left( \begin{array}{l} |
198 |
tim |
2999 |
F_R \\ |
199 |
|
|
\tau _R \\ |
200 |
tim |
2746 |
\end{array} \right) = - \left( {\begin{array}{*{20}c} |
201 |
tim |
2999 |
{\Xi ^{tt} } & {\Xi ^{rt} } \\ |
202 |
|
|
{\Xi ^{tr} } & {\Xi ^{rr} } \\ |
203 |
tim |
2746 |
\end{array}} \right)\left( \begin{array}{l} |
204 |
tim |
2999 |
v \\ |
205 |
|
|
w \\ |
206 |
tim |
2746 |
\end{array} \right) |
207 |
tim |
2999 |
\] |
208 |
|
|
where $F_r$ is the friction force and $\tau _R$ is the friction |
209 |
|
|
torque. |
210 |
tim |
2746 |
|
211 |
tim |
2999 |
\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} |
212 |
tim |
2746 |
|
213 |
tim |
2999 |
For a spherical particle with slip boundary conditions, the |
214 |
|
|
translational and rotational friction constant can be calculated |
215 |
|
|
from Stoke's law, |
216 |
tim |
2746 |
\[ |
217 |
tim |
2999 |
\Xi ^{tt} = \left( {\begin{array}{*{20}c} |
218 |
|
|
{6\pi \eta R} & 0 & 0 \\ |
219 |
|
|
0 & {6\pi \eta R} & 0 \\ |
220 |
|
|
0 & 0 & {6\pi \eta R} \\ |
221 |
|
|
\end{array}} \right) |
222 |
tim |
2746 |
\] |
223 |
tim |
2999 |
and |
224 |
tim |
2746 |
\[ |
225 |
tim |
2999 |
\Xi ^{rr} = \left( {\begin{array}{*{20}c} |
226 |
|
|
{8\pi \eta R^3 } & 0 & 0 \\ |
227 |
|
|
0 & {8\pi \eta R^3 } & 0 \\ |
228 |
|
|
0 & 0 & {8\pi \eta R^3 } \\ |
229 |
|
|
\end{array}} \right) |
230 |
tim |
2746 |
\] |
231 |
tim |
2999 |
where $\eta$ is the viscosity of the solvent and $R$ is the |
232 |
|
|
hydrodynamic radius. |
233 |
|
|
|
234 |
|
|
Other non-spherical shapes, such as cylinders and ellipsoids, are |
235 |
|
|
widely used as references for developing new hydrodynamics theory, |
236 |
|
|
because their properties can be calculated exactly. In 1936, Perrin |
237 |
|
|
extended Stokes's law to general ellipsoids, also called a triaxial |
238 |
|
|
ellipsoid, which is given in Cartesian coordinates |
239 |
|
|
by\cite{Perrin1934, Perrin1936} |
240 |
tim |
2746 |
\[ |
241 |
tim |
2999 |
\frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2 |
242 |
|
|
}} = 1 |
243 |
tim |
2746 |
\] |
244 |
tim |
2999 |
where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately, |
245 |
|
|
due to the complexity of the elliptic integral, only the ellipsoid |
246 |
|
|
with the restriction of two axes being equal, \textit{i.e.} |
247 |
|
|
prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved |
248 |
|
|
exactly. Introducing an elliptic integral parameter $S$ for prolate |
249 |
|
|
ellipsoids : |
250 |
|
|
\[ |
251 |
|
|
S = \frac{2}{{\sqrt {a^2 - b^2 } }}\ln \frac{{a + \sqrt {a^2 - b^2 |
252 |
|
|
} }}{b}, |
253 |
|
|
\] |
254 |
|
|
and oblate ellipsoids: |
255 |
|
|
\[ |
256 |
|
|
S = \frac{2}{{\sqrt {b^2 - a^2 } }}arctg\frac{{\sqrt {b^2 - a^2 } |
257 |
|
|
}}{a}, |
258 |
|
|
\] |
259 |
|
|
one can write down the translational and rotational resistance |
260 |
|
|
tensors |
261 |
|
|
\begin{eqnarray*} |
262 |
|
|
\Xi _a^{tt} & = & 16\pi \eta \frac{{a^2 - b^2 }}{{(2a^2 - b^2 )S - 2a}}. \\ |
263 |
|
|
\Xi _b^{tt} & = & \Xi _c^{tt} = 32\pi \eta \frac{{a^2 - b^2 }}{{(2a^2 - 3b^2 )S + |
264 |
|
|
2a}}, |
265 |
|
|
\end{eqnarray*} |
266 |
|
|
and |
267 |
|
|
\begin{eqnarray*} |
268 |
|
|
\Xi _a^{rr} & = & \frac{{32\pi }}{3}\eta \frac{{(a^2 - b^2 )b^2 }}{{2a - b^2 S}}, \\ |
269 |
|
|
\Xi _b^{rr} & = & \Xi _c^{rr} = \frac{{32\pi }}{3}\eta \frac{{(a^4 - b^4 )}}{{(2a^2 - b^2 )S - 2a}}. |
270 |
|
|
\end{eqnarray*} |
271 |
tim |
2746 |
|
272 |
tim |
2999 |
\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}} |
273 |
|
|
|
274 |
|
|
Unlike spherical and other simply shaped molecules, there is no |
275 |
|
|
analytical solution for the friction tensor for arbitrarily shaped |
276 |
|
|
rigid molecules. The ellipsoid of revolution model and general |
277 |
|
|
triaxial ellipsoid model have been used to approximate the |
278 |
|
|
hydrodynamic properties of rigid bodies. However, since the mapping |
279 |
|
|
from all possible ellipsoidal spaces, $r$-space, to all possible |
280 |
|
|
combination of rotational diffusion coefficients, $D$-space, is not |
281 |
|
|
unique\cite{Wegener1979} as well as the intrinsic coupling between |
282 |
|
|
translational and rotational motion of rigid bodies, general |
283 |
|
|
ellipsoids are not always suitable for modeling arbitrarily shaped |
284 |
|
|
rigid molecules. A number of studies have been devoted to |
285 |
|
|
determining the friction tensor for irregularly shaped rigid bodies |
286 |
|
|
using more advanced methods where the molecule of interest was |
287 |
|
|
modeled by a combinations of spheres\cite{Carrasco1999} and the |
288 |
|
|
hydrodynamics properties of the molecule can be calculated using the |
289 |
|
|
hydrodynamic interaction tensor. Let us consider a rigid assembly of |
290 |
|
|
$N$ beads immersed in a continuous medium. Due to hydrodynamic |
291 |
|
|
interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different |
292 |
|
|
than its unperturbed velocity $v_i$, |
293 |
tim |
2746 |
\[ |
294 |
tim |
2999 |
v'_i = v_i - \sum\limits_{j \ne i} {T_{ij} F_j } |
295 |
tim |
2746 |
\] |
296 |
tim |
2999 |
where $F_i$ is the frictional force, and $T_{ij}$ is the |
297 |
|
|
hydrodynamic interaction tensor. The friction force of $i$th bead is |
298 |
|
|
proportional to its ``net'' velocity |
299 |
tim |
2746 |
\begin{equation} |
300 |
tim |
2999 |
F_i = \zeta _i v_i - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }. |
301 |
|
|
\label{introEquation:tensorExpression} |
302 |
tim |
2746 |
\end{equation} |
303 |
tim |
2999 |
This equation is the basis for deriving the hydrodynamic tensor. In |
304 |
|
|
1930, Oseen and Burgers gave a simple solution to |
305 |
|
|
Eq.~\ref{introEquation:tensorExpression} |
306 |
tim |
2746 |
\begin{equation} |
307 |
tim |
2999 |
T_{ij} = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij} |
308 |
|
|
R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor} |
309 |
tim |
2746 |
\end{equation} |
310 |
tim |
2999 |
Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$. |
311 |
|
|
A second order expression for element of different size was |
312 |
|
|
introduced by Rotne and Prager\cite{Rotne1969} and improved by |
313 |
|
|
Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977} |
314 |
tim |
2746 |
\begin{equation} |
315 |
tim |
2999 |
T_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I + |
316 |
|
|
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma |
317 |
|
|
_i^2 + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} - |
318 |
|
|
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. |
319 |
|
|
\label{introEquation:RPTensorNonOverlapped} |
320 |
tim |
2746 |
\end{equation} |
321 |
tim |
2999 |
Both of the Eq.~\ref{introEquation:oseenTensor} and |
322 |
|
|
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption |
323 |
|
|
$R_{ij} \ge \sigma _i + \sigma _j$. An alternative expression for |
324 |
|
|
overlapping beads with the same radius, $\sigma$, is given by |
325 |
tim |
2746 |
\begin{equation} |
326 |
tim |
2999 |
T_{ij} = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 - |
327 |
|
|
\frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I + |
328 |
|
|
\frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right] |
329 |
|
|
\label{introEquation:RPTensorOverlapped} |
330 |
tim |
2746 |
\end{equation} |
331 |
tim |
2999 |
To calculate the resistance tensor at an arbitrary origin $O$, we |
332 |
|
|
construct a $3N \times 3N$ matrix consisting of $N \times N$ |
333 |
|
|
$B_{ij}$ blocks |
334 |
|
|
\begin{equation} |
335 |
|
|
B = \left( {\begin{array}{*{20}c} |
336 |
|
|
{B_{11} } & \ldots & {B_{1N} } \\ |
337 |
|
|
\vdots & \ddots & \vdots \\ |
338 |
|
|
{B_{N1} } & \cdots & {B_{NN} } \\ |
339 |
|
|
\end{array}} \right), |
340 |
|
|
\end{equation} |
341 |
|
|
where $B_{ij}$ is given by |
342 |
tim |
2746 |
\[ |
343 |
tim |
2999 |
B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij} |
344 |
|
|
)T_{ij} |
345 |
tim |
2746 |
\] |
346 |
tim |
2999 |
where $\delta _{ij}$ is the Kronecker delta function. Inverting the |
347 |
|
|
$B$ matrix, we obtain |
348 |
tim |
2746 |
\[ |
349 |
tim |
2999 |
C = B^{ - 1} = \left( {\begin{array}{*{20}c} |
350 |
|
|
{C_{11} } & \ldots & {C_{1N} } \\ |
351 |
|
|
\vdots & \ddots & \vdots \\ |
352 |
|
|
{C_{N1} } & \cdots & {C_{NN} } \\ |
353 |
|
|
\end{array}} \right), |
354 |
tim |
2746 |
\] |
355 |
tim |
2999 |
which can be partitioned into $N \times N$ $3 \times 3$ block |
356 |
|
|
$C_{ij}$. With the help of $C_{ij}$ and the skew matrix $U_i$ |
357 |
tim |
2746 |
\[ |
358 |
tim |
2999 |
U_i = \left( {\begin{array}{*{20}c} |
359 |
|
|
0 & { - z_i } & {y_i } \\ |
360 |
|
|
{z_i } & 0 & { - x_i } \\ |
361 |
|
|
{ - y_i } & {x_i } & 0 \\ |
362 |
|
|
\end{array}} \right) |
363 |
tim |
2746 |
\] |
364 |
tim |
2999 |
where $x_i$, $y_i$, $z_i$ are the components of the vector joining |
365 |
|
|
bead $i$ and origin $O$, the elements of resistance tensor at |
366 |
|
|
arbitrary origin $O$ can be written as |
367 |
|
|
\begin{eqnarray} |
368 |
|
|
\Xi _{}^{tt} & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\ |
369 |
|
|
\Xi _{}^{tr} & = & \Xi _{}^{rt} = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\ |
370 |
|
|
\Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\ |
371 |
|
|
\label{introEquation:ResistanceTensorArbitraryOrigin} |
372 |
|
|
\end{eqnarray} |
373 |
|
|
The resistance tensor depends on the origin to which they refer. The |
374 |
|
|
proper location for applying the friction force is the center of |
375 |
|
|
resistance (or center of reaction), at which the trace of rotational |
376 |
|
|
resistance tensor, $ \Xi ^{rr}$ reaches a minimum value. |
377 |
|
|
Mathematically, the center of resistance is defined as an unique |
378 |
|
|
point of the rigid body at which the translation-rotation coupling |
379 |
|
|
tensors are symmetric, |
380 |
|
|
\begin{equation} |
381 |
|
|
\Xi^{tr} = \left( {\Xi^{tr} } \right)^T |
382 |
|
|
\label{introEquation:definitionCR} |
383 |
|
|
\end{equation} |
384 |
|
|
From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, |
385 |
|
|
we can easily derive that the translational resistance tensor is |
386 |
|
|
origin independent, while the rotational resistance tensor and |
387 |
|
|
translation-rotation coupling resistance tensor depend on the |
388 |
|
|
origin. Given the resistance tensor at an arbitrary origin $O$, and |
389 |
|
|
a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can |
390 |
|
|
obtain the resistance tensor at $P$ by |
391 |
|
|
\begin{equation} |
392 |
|
|
\begin{array}{l} |
393 |
|
|
\Xi _P^{tt} = \Xi _O^{tt} \\ |
394 |
|
|
\Xi _P^{tr} = \Xi _P^{rt} = \Xi _O^{tr} - U_{OP} \Xi _O^{tt} \\ |
395 |
|
|
\Xi _P^{rr} = \Xi _O^{rr} - U_{OP} \Xi _O^{tt} U_{OP} + \Xi _O^{tr} U_{OP} - U_{OP} \Xi _O^{{tr} ^{^T }} \\ |
396 |
|
|
\end{array} |
397 |
|
|
\label{introEquation:resistanceTensorTransformation} |
398 |
|
|
\end{equation} |
399 |
|
|
where |
400 |
tim |
2746 |
\[ |
401 |
tim |
2999 |
U_{OP} = \left( {\begin{array}{*{20}c} |
402 |
|
|
0 & { - z_{OP} } & {y_{OP} } \\ |
403 |
|
|
{z_i } & 0 & { - x_{OP} } \\ |
404 |
|
|
{ - y_{OP} } & {x_{OP} } & 0 \\ |
405 |
|
|
\end{array}} \right) |
406 |
tim |
2746 |
\] |
407 |
tim |
2999 |
Using Eq.~\ref{introEquation:definitionCR} and |
408 |
|
|
Eq.~\ref{introEquation:resistanceTensorTransformation}, one can |
409 |
|
|
locate the position of center of resistance, |
410 |
|
|
\begin{eqnarray*} |
411 |
|
|
\left( \begin{array}{l} |
412 |
|
|
x_{OR} \\ |
413 |
|
|
y_{OR} \\ |
414 |
|
|
z_{OR} \\ |
415 |
|
|
\end{array} \right) & = &\left( {\begin{array}{*{20}c} |
416 |
|
|
{(\Xi _O^{rr} )_{yy} + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} } \\ |
417 |
|
|
{ - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz} + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} } \\ |
418 |
|
|
{ - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx} + (\Xi _O^{rr} )_{yy} } \\ |
419 |
|
|
\end{array}} \right)^{ - 1} \\ |
420 |
|
|
& & \left( \begin{array}{l} |
421 |
|
|
(\Xi _O^{tr} )_{yz} - (\Xi _O^{tr} )_{zy} \\ |
422 |
|
|
(\Xi _O^{tr} )_{zx} - (\Xi _O^{tr} )_{xz} \\ |
423 |
|
|
(\Xi _O^{tr} )_{xy} - (\Xi _O^{tr} )_{yx} \\ |
424 |
|
|
\end{array} \right) \\ |
425 |
|
|
\end{eqnarray*} |
426 |
|
|
where $x_OR$, $y_OR$, $z_OR$ are the components of the vector |
427 |
|
|
joining center of resistance $R$ and origin $O$. |
428 |
tim |
2746 |
|
429 |
tim |
2999 |
\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
430 |
tim |
2746 |
|
431 |
tim |
2999 |
Consider the Langevin equations of motion in generalized coordinates |
432 |
tim |
2746 |
\begin{equation} |
433 |
|
|
M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)} + F_{r,i} (t) |
434 |
|
|
\label{LDGeneralizedForm} |
435 |
|
|
\end{equation} |
436 |
|
|
where $M_i$ is a $6\times6$ generalized diagonal mass (include mass |
437 |
|
|
and moment of inertial) matrix and $V_i$ is a generalized velocity, |
438 |
tim |
2999 |
$V_i = V_i(v_i,\omega _i)$. The right side of |
439 |
|
|
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in |
440 |
tim |
2746 |
lab-fixed frame, systematic force $F_{s,i}$, dissipative force |
441 |
|
|
$F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the |
442 |
|
|
system in Newtownian mechanics typically refers to lab-fixed frame, |
443 |
|
|
it is also convenient to handle the rotation of rigid body in |
444 |
|
|
body-fixed frame. Thus the friction and random forces are calculated |
445 |
|
|
in body-fixed frame and converted back to lab-fixed frame by: |
446 |
|
|
\[ |
447 |
|
|
\begin{array}{l} |
448 |
tim |
2999 |
F_{f,i}^l (t) = Q^T F_{f,i}^b (t), \\ |
449 |
|
|
F_{r,i}^l (t) = Q^T F_{r,i}^b (t). \\ |
450 |
|
|
\end{array} |
451 |
tim |
2746 |
\] |
452 |
|
|
Here, the body-fixed friction force $F_{r,i}^b$ is proportional to |
453 |
|
|
the body-fixed velocity at center of resistance $v_{R,i}^b$ and |
454 |
tim |
2999 |
angular velocity $\omega _i$ |
455 |
tim |
2746 |
\begin{equation} |
456 |
|
|
F_{r,i}^b (t) = \left( \begin{array}{l} |
457 |
|
|
f_{r,i}^b (t) \\ |
458 |
|
|
\tau _{r,i}^b (t) \\ |
459 |
|
|
\end{array} \right) = - \left( {\begin{array}{*{20}c} |
460 |
|
|
{\Xi _{R,t} } & {\Xi _{R,c}^T } \\ |
461 |
|
|
{\Xi _{R,c} } & {\Xi _{R,r} } \\ |
462 |
|
|
\end{array}} \right)\left( \begin{array}{l} |
463 |
|
|
v_{R,i}^b (t) \\ |
464 |
|
|
\omega _i (t) \\ |
465 |
|
|
\end{array} \right), |
466 |
|
|
\end{equation} |
467 |
|
|
while the random force $F_{r,i}^l$ is a Gaussian stochastic variable |
468 |
|
|
with zero mean and variance |
469 |
|
|
\begin{equation} |
470 |
|
|
\left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle = |
471 |
|
|
\left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle = |
472 |
tim |
2999 |
2k_B T\Xi _R \delta (t - t'). \label{randomForce} |
473 |
tim |
2746 |
\end{equation} |
474 |
|
|
The equation of motion for $v_i$ can be written as |
475 |
|
|
\begin{equation} |
476 |
|
|
m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) + |
477 |
|
|
f_{r,i}^l (t) |
478 |
|
|
\end{equation} |
479 |
|
|
Since the frictional force is applied at the center of resistance |
480 |
|
|
which generally does not coincide with the center of mass, an extra |
481 |
|
|
torque is exerted at the center of mass. Thus, the net body-fixed |
482 |
|
|
frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is |
483 |
|
|
given by |
484 |
|
|
\begin{equation} |
485 |
|
|
\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b |
486 |
|
|
\end{equation} |
487 |
|
|
where $r_{MR}$ is the vector from the center of mass to the center |
488 |
tim |
2999 |
of the resistance. Instead of integrating the angular velocity in |
489 |
|
|
lab-fixed frame, we consider the equation of angular momentum in |
490 |
|
|
body-fixed frame |
491 |
tim |
2746 |
\begin{equation} |
492 |
tim |
2999 |
\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t) |
493 |
|
|
+ \tau _{r,i}^b(t) |
494 |
tim |
2746 |
\end{equation} |
495 |
|
|
Embedding the friction terms into force and torque, one can |
496 |
|
|
integrate the langevin equations of motion for rigid body of |
497 |
|
|
arbitrary shape in a velocity-Verlet style 2-part algorithm, where |
498 |
|
|
$h= \delta t$: |
499 |
|
|
|
500 |
tim |
2999 |
{\tt moveA:} |
501 |
tim |
2746 |
\begin{align*} |
502 |
tim |
2999 |
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
503 |
|
|
+ \frac{h}{2} \left( {\bf f}(t) / m \right), \\ |
504 |
|
|
% |
505 |
|
|
{\bf r}(t + h) &\leftarrow {\bf r}(t) |
506 |
|
|
+ h {\bf v}\left(t + h / 2 \right), \\ |
507 |
|
|
% |
508 |
|
|
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
509 |
|
|
+ \frac{h}{2} {\bf \tau}^b(t), \\ |
510 |
|
|
% |
511 |
|
|
\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
512 |
|
|
(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right). |
513 |
tim |
2746 |
\end{align*} |
514 |
|
|
In this context, the $\mathrm{rotate}$ function is the reversible |
515 |
tim |
2999 |
product of the three body-fixed rotations, |
516 |
tim |
2746 |
\begin{equation} |
517 |
|
|
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
518 |
|
|
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y |
519 |
|
|
/ 2) \cdot \mathsf{G}_x(a_x /2), |
520 |
|
|
\end{equation} |
521 |
|
|
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, |
522 |
tim |
2999 |
rotates both the rotation matrix ($\mathsf{Q}$) and the body-fixed |
523 |
|
|
angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed |
524 |
|
|
axis $\alpha$, |
525 |
tim |
2746 |
\begin{equation} |
526 |
|
|
\mathsf{G}_\alpha( \theta ) = \left\{ |
527 |
|
|
\begin{array}{lcl} |
528 |
tim |
2999 |
\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
529 |
tim |
2746 |
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf |
530 |
|
|
j}(0). |
531 |
|
|
\end{array} |
532 |
|
|
\right. |
533 |
|
|
\end{equation} |
534 |
|
|
$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis |
535 |
|
|
rotation matrix. For example, in the small-angle limit, the |
536 |
|
|
rotation matrix around the body-fixed x-axis can be approximated as |
537 |
|
|
\begin{equation} |
538 |
|
|
\mathsf{R}_x(\theta) \approx \left( |
539 |
|
|
\begin{array}{ccc} |
540 |
|
|
1 & 0 & 0 \\ |
541 |
|
|
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+ |
542 |
|
|
\theta^2 / 4} \\ |
543 |
|
|
0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + |
544 |
|
|
\theta^2 / 4} |
545 |
|
|
\end{array} |
546 |
|
|
\right). |
547 |
|
|
\end{equation} |
548 |
tim |
2999 |
All other rotations follow in a straightforward manner. After the |
549 |
|
|
first part of the propagation, the forces and body-fixed torques are |
550 |
|
|
calculated at the new positions and orientations |
551 |
tim |
2746 |
|
552 |
tim |
2999 |
{\tt doForces:} |
553 |
|
|
\begin{align*} |
554 |
|
|
{\bf f}(t + h) &\leftarrow |
555 |
|
|
- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\ |
556 |
|
|
% |
557 |
|
|
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) |
558 |
|
|
\times \frac{\partial V}{\partial {\bf u}}, \\ |
559 |
|
|
% |
560 |
|
|
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h) |
561 |
|
|
\cdot {\bf \tau}^s(t + h). |
562 |
|
|
\end{align*} |
563 |
tim |
2746 |
Once the forces and torques have been obtained at the new time step, |
564 |
|
|
the velocities can be advanced to the same time value. |
565 |
|
|
|
566 |
tim |
2999 |
{\tt moveB:} |
567 |
tim |
2746 |
\begin{align*} |
568 |
tim |
2999 |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2 |
569 |
|
|
\right) |
570 |
|
|
+ \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\ |
571 |
|
|
% |
572 |
|
|
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2 |
573 |
|
|
\right) |
574 |
|
|
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
575 |
tim |
2746 |
\end{align*} |
576 |
|
|
|
577 |
gezelter |
3305 |
\section{Results} |
578 |
gezelter |
3302 |
In order to validate our Langevin integrator for arbitrarily-shaped |
579 |
gezelter |
3305 |
rigid bodies, we implemented the algorithm in {\sc |
580 |
|
|
oopse}\cite{Meineke2005} and compared the results of this algorithm |
581 |
|
|
with the known |
582 |
gezelter |
3302 |
hydrodynamic limiting behavior for a few model systems, and to |
583 |
|
|
microcanonical molecular dynamics simulations for some more |
584 |
|
|
complicated bodies. The model systems and their analytical behavior |
585 |
|
|
(if known) are summarized below. Parameters for the primary particles |
586 |
|
|
comprising our model systems are given in table \ref{tab:parameters}, |
587 |
|
|
and a sketch of the arrangement of these primary particles into the |
588 |
gezelter |
3305 |
model rigid bodies is shown in figure \ref{fig:models}. In table |
589 |
|
|
\ref{tab:parameters}, $d$ and $l$ are the physical dimensions of |
590 |
|
|
ellipsoidal (Gay-Berne) particles. For spherical particles, the value |
591 |
|
|
of the Lennard-Jones $\sigma$ parameter is the particle diameter |
592 |
|
|
($d$). Gay-Berne ellipsoids have an energy scaling parameter, |
593 |
|
|
$\epsilon^s$, which describes the well depth for two identical |
594 |
|
|
ellipsoids in a {\it side-by-side} configuration. Additionally, a |
595 |
|
|
well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, |
596 |
|
|
describes the ratio between the well depths in the {\it end-to-end} |
597 |
|
|
and side-by-side configurations. For spheres, $\epsilon^r \equiv 1$. |
598 |
|
|
Moments of inertia are also required to describe the motion of primary |
599 |
|
|
particles with orientational degrees of freedom. |
600 |
gezelter |
3299 |
|
601 |
gezelter |
3302 |
\begin{table*} |
602 |
|
|
\begin{minipage}{\linewidth} |
603 |
|
|
\begin{center} |
604 |
|
|
\caption{Parameters for the primary particles in use by the rigid body |
605 |
|
|
models in figure \ref{fig:models}.} |
606 |
|
|
\begin{tabular}{lrcccccccc} |
607 |
|
|
\hline |
608 |
|
|
& & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\ |
609 |
|
|
& & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ (kcal/mol) & $\epsilon^r$ & |
610 |
|
|
$m$ (amu) & $I_{xx}$ & $I_{yy}$ & $I_{zz}$ \\ \hline |
611 |
gezelter |
3308 |
Sphere & & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
612 |
gezelter |
3302 |
Ellipsoid & & 4.6 & 13.8 & 0.8 & 0.2 & 200 & 2105 & 2105 & 421 \\ |
613 |
gezelter |
3308 |
Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
614 |
gezelter |
3302 |
Banana &(3 identical ellipsoids)& 4.2 & 11.2 & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\ |
615 |
|
|
Lipid: & Spherical Head & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\ |
616 |
|
|
& Ellipsoidal Tail & 4.6 & 13.8 & 0.8 & 0.2 & 760 & 45000 & 45000 & 9000 \\ |
617 |
|
|
Solvent & & 4.7 & $= d$ & 0.8 & 1 & 72.06 & & & \\ |
618 |
|
|
\hline |
619 |
|
|
\end{tabular} |
620 |
|
|
\label{tab:parameters} |
621 |
|
|
\end{center} |
622 |
|
|
\end{minipage} |
623 |
|
|
\end{table*} |
624 |
|
|
|
625 |
gezelter |
3305 |
\begin{figure} |
626 |
|
|
\centering |
627 |
|
|
\includegraphics[width=3in]{sketch} |
628 |
|
|
\caption[Sketch of the model systems]{A sketch of the model systems |
629 |
|
|
used in evaluating the behavior of the rigid body Langevin |
630 |
|
|
integrator.} \label{fig:models} |
631 |
|
|
\end{figure} |
632 |
|
|
|
633 |
gezelter |
3302 |
\subsection{Simulation Methodology} |
634 |
|
|
|
635 |
|
|
We performed reference microcanonical simulations with explicit |
636 |
|
|
solvents for each of the different model system. In each case there |
637 |
|
|
was one solute model and 1929 solvent molecules present in the |
638 |
|
|
simulation box. All simulations were equilibrated using a |
639 |
|
|
constant-pressure and temperature integrator with target values of 300 |
640 |
|
|
K for the temperature and 1 atm for pressure. Following this stage, |
641 |
|
|
further equilibration and sampling was done in a microcanonical |
642 |
gezelter |
3305 |
ensemble. Since the model bodies are typically quite massive, we were |
643 |
|
|
able to use a time step of 25 fs. A switching function was applied to |
644 |
gezelter |
3302 |
all potentials to smoothly turn off the interactions between a range |
645 |
|
|
of $22$ and $25$ \AA. The switching function was the standard (cubic) |
646 |
|
|
function, |
647 |
|
|
\begin{equation} |
648 |
|
|
s(r) = |
649 |
|
|
\begin{cases} |
650 |
|
|
1 & \text{if $r \le r_{\text{sw}}$},\\ |
651 |
|
|
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
652 |
|
|
{(r_{\text{cut}} - r_{\text{sw}})^3} |
653 |
|
|
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
654 |
|
|
0 & \text{if $r > r_{\text{cut}}$.} |
655 |
|
|
\end{cases} |
656 |
|
|
\label{eq:switchingFunc} |
657 |
|
|
\end{equation} |
658 |
|
|
To measure shear viscosities from our microcanonical simulations, we |
659 |
|
|
used the Einstein form of the pressure correlation function,\cite{hess:209} |
660 |
|
|
\begin{equation} |
661 |
|
|
\eta = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \langle \left( |
662 |
|
|
\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \rangle_{t_0}. |
663 |
|
|
\label{eq:shear} |
664 |
|
|
\end{equation} |
665 |
|
|
A similar form exists for the bulk viscosity |
666 |
|
|
\begin{equation} |
667 |
|
|
\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \langle \left( |
668 |
|
|
\int_{t_0}^{t_0 + t} |
669 |
|
|
\left(P\left(t'\right)-\langle P \rangle \right)dt' |
670 |
|
|
\right)^2 \rangle_{t_0}. |
671 |
|
|
\end{equation} |
672 |
|
|
Alternatively, the shear viscosity can also be calculated using a |
673 |
|
|
Green-Kubo formula with the off-diagonal pressure tensor correlation function, |
674 |
|
|
\begin{equation} |
675 |
|
|
\eta = \frac{V}{k_B T} \int_0^{\infty} \langle P_{xz}(t_0) P_{xz}(t_0 |
676 |
|
|
+ t) \rangle_{t_0} dt, |
677 |
|
|
\end{equation} |
678 |
|
|
although this method converges extremely slowly and is not practical |
679 |
|
|
for obtaining viscosities from molecular dynamics simulations. |
680 |
|
|
|
681 |
|
|
The Langevin dynamics for the different model systems were performed |
682 |
|
|
at the same temperature as the average temperature of the |
683 |
|
|
microcanonical simulations and with a solvent viscosity taken from |
684 |
gezelter |
3305 |
Eq. (\ref{eq:shear}) applied to these simulations. We used 1024 |
685 |
|
|
independent solute simulations to obtain statistics on our Langevin |
686 |
|
|
integrator. |
687 |
gezelter |
3302 |
|
688 |
|
|
\subsection{Analysis} |
689 |
|
|
|
690 |
|
|
The quantities of interest when comparing the Langevin integrator to |
691 |
|
|
analytic hydrodynamic equations and to molecular dynamics simulations |
692 |
|
|
are typically translational diffusion constants and orientational |
693 |
|
|
relaxation times. Translational diffusion constants for point |
694 |
|
|
particles are computed easily from the long-time slope of the |
695 |
|
|
mean-square displacement, |
696 |
|
|
\begin{equation} |
697 |
|
|
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle, |
698 |
|
|
\end{equation} |
699 |
|
|
of the solute molecules. For models in which the translational |
700 |
gezelter |
3305 |
diffusion tensor (${\bf D}_{tt}$) has non-degenerate eigenvalues |
701 |
|
|
(i.e. any non-spherically-symmetric rigid body), it is possible to |
702 |
|
|
compute the diffusive behavior for motion parallel to each body-fixed |
703 |
|
|
axis by projecting the displacement of the particle onto the |
704 |
|
|
body-fixed reference frame at $t=0$. With an isotropic solvent, as we |
705 |
|
|
have used in this study, there are differences between the three |
706 |
gezelter |
3302 |
diffusion constants, but these must converge to the same value at |
707 |
|
|
longer times. Translational diffusion constants for the different |
708 |
gezelter |
3305 |
shaped models are shown in table \ref{tab:translation}. |
709 |
gezelter |
3302 |
|
710 |
gezelter |
3305 |
In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational |
711 |
gezelter |
3302 |
diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object |
712 |
|
|
{\it around} a particular body-fixed axis and {\it not} the diffusion |
713 |
|
|
of a vector pointing along the axis. However, these eigenvalues can |
714 |
|
|
be combined to find 5 characteristic rotational relaxation |
715 |
gezelter |
3305 |
times,\cite{PhysRev.119.53,Berne90} |
716 |
gezelter |
3302 |
\begin{eqnarray} |
717 |
gezelter |
3305 |
1 / \tau_1 & = & 6 D_r + 2 \Delta \\ |
718 |
|
|
1 / \tau_2 & = & 6 D_r - 2 \Delta \\ |
719 |
|
|
1 / \tau_3 & = & 3 (D_r + D_1) \\ |
720 |
|
|
1 / \tau_4 & = & 3 (D_r + D_2) \\ |
721 |
|
|
1 / \tau_5 & = & 3 (D_r + D_3) |
722 |
gezelter |
3302 |
\end{eqnarray} |
723 |
|
|
where |
724 |
|
|
\begin{equation} |
725 |
|
|
D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right) |
726 |
|
|
\end{equation} |
727 |
|
|
and |
728 |
|
|
\begin{equation} |
729 |
gezelter |
3305 |
\Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2} |
730 |
gezelter |
3302 |
\end{equation} |
731 |
gezelter |
3305 |
Each of these characteristic times can be used to predict the decay of |
732 |
|
|
part of the rotational correlation function when $\ell = 2$, |
733 |
gezelter |
3302 |
\begin{equation} |
734 |
gezelter |
3305 |
C_2(t) = \frac{a^2}{N^2} e^{-t/\tau_1} + \frac{b^2}{N^2} e^{-t/\tau_2}. |
735 |
gezelter |
3302 |
\end{equation} |
736 |
gezelter |
3305 |
This is the same as the $F^2_{0,0}(t)$ correlation function that |
737 |
|
|
appears in Ref. \citen{Berne90}. The amplitudes of the two decay |
738 |
|
|
terms are expressed in terms of three dimensionless functions of the |
739 |
|
|
eigenvalues: $a = \sqrt{3} (D_1 - D_2)$, $b = (2D_3 - D_1 - D_2 + |
740 |
|
|
2\Delta)$, and $N = 2 \sqrt{\Delta b}$. Similar expressions can be |
741 |
|
|
obtained for other angular momentum correlation |
742 |
|
|
functions.\cite{PhysRev.119.53,Berne90} In all of the model systems we |
743 |
|
|
studied, only one of the amplitudes of the two decay terms was |
744 |
|
|
non-zero, so it was possible to derive a single relaxation time for |
745 |
|
|
each of the hydrodynamic tensors. In many cases, these characteristic |
746 |
|
|
times are averaged and reported in the literature as a single relaxation |
747 |
|
|
time,\cite{Garcia-de-la-Torre:1997qy} |
748 |
gezelter |
3302 |
\begin{equation} |
749 |
gezelter |
3305 |
1 / \tau_0 = \frac{1}{5} \sum_{i=1}^5 \tau_{i}^{-1}, |
750 |
|
|
\end{equation} |
751 |
|
|
although for the cases reported here, this averaging is not necessary |
752 |
|
|
and only one of the five relaxation times is relevant. |
753 |
|
|
|
754 |
|
|
To test the Langevin integrator's behavior for rotational relaxation, |
755 |
|
|
we have compared the analytical orientational relaxation times (if |
756 |
|
|
they are known) with the general result from the diffusion tensor and |
757 |
|
|
with the results from both the explicitly solvated molecular dynamics |
758 |
|
|
and Langevin simulations. Relaxation times from simulations (both |
759 |
|
|
microcanonical and Langevin), were computed using Legendre polynomial |
760 |
|
|
correlation functions for a unit vector (${\bf u}$) fixed along one or |
761 |
|
|
more of the body-fixed axes of the model. |
762 |
|
|
\begin{equation} |
763 |
gezelter |
3302 |
C_{\ell}(t) = \langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf |
764 |
|
|
u}_{i}(0) \right) |
765 |
|
|
\end{equation} |
766 |
|
|
For simulations in the high-friction limit, orientational correlation |
767 |
|
|
times can then be obtained from exponential fits of this function, or by |
768 |
|
|
integrating, |
769 |
|
|
\begin{equation} |
770 |
gezelter |
3305 |
\tau = \ell (\ell + 1) \int_0^{\infty} C_{\ell}(t) dt. |
771 |
gezelter |
3302 |
\end{equation} |
772 |
gezelter |
3305 |
In lower-friction solvents, the Legendre correlation functions often |
773 |
|
|
exhibit non-exponential decay, and may not be characterized by a |
774 |
|
|
single decay constant. |
775 |
gezelter |
3302 |
|
776 |
|
|
In table \ref{tab:rotation} we show the characteristic rotational |
777 |
|
|
relaxation times (based on the diffusion tensor) for each of the model |
778 |
|
|
systems compared with the values obtained via microcanonical and Langevin |
779 |
|
|
simulations. |
780 |
|
|
|
781 |
gezelter |
3305 |
\subsection{Spherical particles} |
782 |
gezelter |
3299 |
Our model system for spherical particles was a Lennard-Jones sphere of |
783 |
|
|
diameter ($\sigma$) 6.5 \AA\ in a sea of smaller spheres ($\sigma$ = |
784 |
|
|
4.7 \AA). The well depth ($\epsilon$) for both particles was set to |
785 |
gezelter |
3302 |
an arbitrary value of 0.8 kcal/mol. |
786 |
gezelter |
3299 |
|
787 |
|
|
The Stokes-Einstein behavior of large spherical particles in |
788 |
|
|
hydrodynamic flows is well known, giving translational friction |
789 |
|
|
coefficients of $6 \pi \eta R$ (stick boundary conditions) and |
790 |
gezelter |
3302 |
rotational friction coefficients of $8 \pi \eta R^3$. Recently, |
791 |
|
|
Schmidt and Skinner have computed the behavior of spherical tag |
792 |
|
|
particles in molecular dynamics simulations, and have shown that {\it |
793 |
|
|
slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more |
794 |
gezelter |
3299 |
appropriate for molecule-sized spheres embedded in a sea of spherical |
795 |
gezelter |
3305 |
qsolvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
796 |
gezelter |
3299 |
|
797 |
|
|
Our simulation results show similar behavior to the behavior observed |
798 |
gezelter |
3302 |
by Schmidt and Skinner. The diffusion constant obtained from our |
799 |
gezelter |
3299 |
microcanonical molecular dynamics simulations lies between the slip |
800 |
|
|
and stick boundary condition results obtained via Stokes-Einstein |
801 |
|
|
behavior. Since the Langevin integrator assumes Stokes-Einstein stick |
802 |
|
|
boundary conditions in calculating the drag and random forces for |
803 |
|
|
spherical particles, our Langevin routine obtains nearly quantitative |
804 |
|
|
agreement with the hydrodynamic results for spherical particles. One |
805 |
|
|
avenue for improvement of the method would be to compute elements of |
806 |
|
|
$\Xi_{tt}$ assuming behavior intermediate between the two boundary |
807 |
gezelter |
3302 |
conditions. |
808 |
gezelter |
3299 |
|
809 |
|
|
In these simulations, our spherical particles were structureless, so |
810 |
|
|
there is no way to obtain rotational correlation times for this model |
811 |
|
|
system. |
812 |
|
|
|
813 |
|
|
\subsubsection{Ellipsoids} |
814 |
|
|
For uniaxial ellipsoids ($a > b = c$) , Perrin's formulae for both |
815 |
|
|
translational and rotational diffusion of each of the body-fixed axes |
816 |
|
|
can be combined to give a single translational diffusion |
817 |
gezelter |
3302 |
constant,\cite{Berne90} |
818 |
gezelter |
3299 |
\begin{equation} |
819 |
|
|
D = \frac{k_B T}{6 \pi \eta a} G(\rho), |
820 |
|
|
\label{Dperrin} |
821 |
|
|
\end{equation} |
822 |
|
|
as well as a single rotational diffusion coefficient, |
823 |
|
|
\begin{equation} |
824 |
|
|
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2) |
825 |
|
|
G(\rho) - 1}{1 - \rho^4} \right\}. |
826 |
|
|
\label{ThetaPerrin} |
827 |
|
|
\end{equation} |
828 |
|
|
In these expressions, $G(\rho)$ is a function of the axial ratio |
829 |
|
|
($\rho = b / a$), which for prolate ellipsoids, is |
830 |
|
|
\begin{equation} |
831 |
|
|
G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 - |
832 |
|
|
\rho^2)^{1/2}}{\rho} \right\} |
833 |
|
|
\label{GPerrin} |
834 |
|
|
\end{equation} |
835 |
|
|
Again, there is some uncertainty about the correct boundary conditions |
836 |
|
|
to use for molecular-scale ellipsoids in a sea of similarly-sized |
837 |
|
|
solvent particles. Ravichandran and Bagchi found that {\it slip} |
838 |
gezelter |
3302 |
boundary conditions most closely resembled the simulation |
839 |
|
|
results,\cite{Ravichandran:1999fk} in agreement with earlier work of |
840 |
|
|
Tang and Evans.\cite{TANG:1993lr} |
841 |
gezelter |
3299 |
|
842 |
gezelter |
3305 |
Even though there are analytic resistance tensors for ellipsoids, we |
843 |
|
|
constructed a rough-shell model using 2135 beads (each with a diameter |
844 |
|
|
of 0.25 \AA) to approximate the shape of the modle ellipsoid. We |
845 |
|
|
compared the Langevin dynamics from both the simple ellipsoidal |
846 |
|
|
resistance tensor and the rough shell approximation with |
847 |
|
|
microcanonical simulations and the predictions of Perrin. As in the |
848 |
|
|
case of our spherical model system, the Langevin integrator reproduces |
849 |
|
|
almost exactly the behavior of the Perrin formulae (which is |
850 |
|
|
unsurprising given that the Perrin formulae were used to derive the |
851 |
gezelter |
3299 |
drag and random forces applied to the ellipsoid). We obtain |
852 |
|
|
translational diffusion constants and rotational correlation times |
853 |
|
|
that are within a few percent of the analytic values for both the |
854 |
|
|
exact treatment of the diffusion tensor as well as the rough-shell |
855 |
|
|
model for the ellipsoid. |
856 |
|
|
|
857 |
gezelter |
3308 |
The translational diffusion constants from the microcanonical simulations |
858 |
|
|
agree well with the predictions of the Perrin model, although the rotational |
859 |
|
|
correlation times are a factor of 2 shorter than expected from hydrodynamic |
860 |
|
|
theory. One explanation for the slower rotation |
861 |
|
|
of explicitly-solvated ellipsoids is the possibility that solute-solvent |
862 |
|
|
collisions happen at both ends of the solute whenever the principal |
863 |
|
|
axis of the ellipsoid is turning. In the upper portion of figure |
864 |
|
|
\ref{fig:explanation} we sketch a physical picture of this explanation. |
865 |
|
|
Since our Langevin integrator is providing nearly quantitative agreement with |
866 |
|
|
the Perrin model, it also predicts orientational diffusion for ellipsoids that |
867 |
|
|
exceed explicitly solvated correlation times by a factor of two. |
868 |
gezelter |
3299 |
|
869 |
gezelter |
3308 |
\begin{figure} |
870 |
|
|
\centering |
871 |
|
|
\includegraphics[width=6in]{explanation} |
872 |
|
|
\caption[Explanations of the differences between orientational correlation times for explicitly-solvated models and hydrodynamics predictions]{Explanations of the differences between orientational correlation times for explicitly-solvated models and hydrodynamic predictions. For the ellipsoids (upper figures), rotation of the principal axis can involve correlated collisions at both sides of the solute. In the rigid dumbbell model (lower figures), the large size of the explicit solvent spheres prevents them from coming in contact with a substantial fraction of the surface area of the dumbbell. Therefore, the explicit solvent only provides drag over a substantially reduced surface area of this model, where the hydrodynamic theories utilize the entire surface area for estimating rotational diffusion. |
873 |
|
|
} \label{fig:explanation} |
874 |
|
|
\end{figure} |
875 |
|
|
|
876 |
gezelter |
3302 |
\subsubsection{Rigid dumbbells} |
877 |
|
|
Perhaps the only {\it composite} rigid body for which analytic |
878 |
|
|
expressions for the hydrodynamic tensor are available is the |
879 |
|
|
two-sphere dumbbell model. This model consists of two non-overlapping |
880 |
|
|
spheres held by a rigid bond connecting their centers. There are |
881 |
|
|
competing expressions for the 6x6 resistance tensor for this |
882 |
|
|
model. Equation (\ref{introEquation:oseenTensor}) above gives the |
883 |
|
|
original Oseen tensor, while the second order expression introduced by |
884 |
|
|
Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la |
885 |
|
|
Torre and Bloomfield,\cite{Torre1977} is given above as |
886 |
gezelter |
3299 |
Eq. (\ref{introEquation:RPTensorNonOverlapped}). In our case, we use |
887 |
|
|
a model dumbbell in which the two spheres are identical Lennard-Jones |
888 |
|
|
particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at |
889 |
gezelter |
3302 |
a distance of 6.532 \AA. |
890 |
gezelter |
3299 |
|
891 |
|
|
The theoretical values for the translational diffusion constant of the |
892 |
|
|
dumbbell are calculated from the work of Stimson and Jeffery, who |
893 |
|
|
studied the motion of this system in a flow parallel to the |
894 |
gezelter |
3302 |
inter-sphere axis,\cite{Stimson:1926qy} and Davis, who studied the |
895 |
|
|
motion in a flow {\it perpendicular} to the inter-sphere |
896 |
|
|
axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it |
897 |
|
|
orientational} correlation times for this model system (other than |
898 |
gezelter |
3305 |
those derived from the 6 x 6 tensors mentioned above). |
899 |
gezelter |
3299 |
|
900 |
gezelter |
3305 |
The bead model for this model system comprises the two large spheres |
901 |
|
|
by themselves, while the rough shell approximation used 3368 separate |
902 |
|
|
beads (each with a diameter of 0.25 \AA) to approximate the shape of |
903 |
|
|
the rigid body. The hydrodynamics tensors computed from both the bead |
904 |
|
|
and rough shell models are remarkably similar. Computing the initial |
905 |
|
|
hydrodynamic tensor for a rough shell model can be quite expensive (in |
906 |
|
|
this case it requires inverting a 10104 x 10104 matrix), while the |
907 |
|
|
bead model is typically easy to compute (in this case requiring |
908 |
gezelter |
3308 |
inversion of a 6 x 6 matrix). |
909 |
gezelter |
3305 |
|
910 |
gezelter |
3308 |
\begin{figure} |
911 |
|
|
\centering |
912 |
|
|
\includegraphics[width=3in]{RoughShell} |
913 |
|
|
\caption[Model rigid bodies and their rough shell approximations]{The |
914 |
|
|
model rigid bodies (left column) used to test this algorithm and their |
915 |
|
|
rough-shell approximations (right-column) that were used to compute |
916 |
|
|
the hydrodynamic tensors. The top two models (ellipsoid and dumbbell) |
917 |
|
|
have analytic solutions and were used to test the rough shell |
918 |
|
|
approximation. The lower two models (banana and lipid) were compared |
919 |
|
|
with explicitly-solvated molecular dynamics simulations. } |
920 |
|
|
\label{fig:roughShell} |
921 |
|
|
\end{figure} |
922 |
|
|
|
923 |
|
|
|
924 |
gezelter |
3305 |
Once the hydrodynamic tensor has been computed, there is no additional |
925 |
|
|
penalty for carrying out a Langevin simulation with either of the two |
926 |
|
|
different hydrodynamics models. Our naive expectation is that since |
927 |
|
|
the rigid body's surface is roughened under the various shell models, |
928 |
|
|
the diffusion constants will be even farther from the ``slip'' |
929 |
|
|
boundary conditions than observed for the bead model (which uses a |
930 |
|
|
Stokes-Einstein model to arrive at the hydrodynamic tensor). For the |
931 |
|
|
dumbbell, this prediction is correct although all of the Langevin |
932 |
|
|
diffusion constants are within 6\% of the diffusion constant predicted |
933 |
|
|
from the fully solvated system. |
934 |
|
|
|
935 |
gezelter |
3308 |
For rotational motion, Langevin integration (and the hydrodynamic tensor) |
936 |
|
|
yields rotational correlation times that are substantially shorter than those |
937 |
|
|
obtained from explicitly-solvated simulations. It is likely that this is due |
938 |
|
|
to the large size of the explicit solvent spheres, a feature that prevents |
939 |
|
|
the solvent from coming in contact with a substantial fraction of the surface |
940 |
|
|
area of the dumbbell. Therefore, the explicit solvent only provides drag |
941 |
|
|
over a substantially reduced surface area of this model, while the |
942 |
|
|
hydrodynamic theories utilize the entire surface area for estimating |
943 |
|
|
rotational diffusion. A sketch of the free volume available in the explicit |
944 |
|
|
solvent simulations is shown in figure \ref{fig:explanation}. |
945 |
gezelter |
3305 |
|
946 |
gezelter |
3299 |
\subsubsection{Ellipsoidal-composite banana-shaped molecules} |
947 |
gezelter |
3308 |
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids |
948 |
|
|
have been used by Orlandi {\it et al.} to observe |
949 |
gezelter |
3299 |
mesophases in coarse-grained models bent-core liquid crystalline |
950 |
gezelter |
3308 |
molecules.\cite{Orlandi:2006fk} We have used the same overlapping |
951 |
gezelter |
3299 |
ellipsoids as a way to test the behavior of our algorithm for a |
952 |
|
|
structure of some interest to the materials science community, |
953 |
|
|
although since we are interested in capturing only the hydrodynamic |
954 |
gezelter |
3308 |
behavior of this model, we have left out the dipolar interactions of the |
955 |
gezelter |
3299 |
original Orlandi model. |
956 |
gezelter |
3308 |
|
957 |
|
|
A reference system composed of a single banana rigid body embedded in a |
958 |
|
|
sea of 1929 solvent particles was created and run under standard |
959 |
|
|
(microcanonical) molecular dynamics. The resulting viscosity of this |
960 |
|
|
mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})). |
961 |
|
|
To calculate the hydrodynamic properties of the banana rigid body model, |
962 |
|
|
we created a rough shell (see Fig.~\ref{langevin:roughShell}), in which |
963 |
|
|
the banana is represented as a ``shell'' made of 3321 identical beads |
964 |
|
|
(0.25 \AA in diameter) distributed on the surface. Applying the |
965 |
|
|
procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
966 |
|
|
identified the center of resistance at (0 $\rm{\AA}$, 0.81 $\rm{\AA}$, 0 $\rm{\AA}$), as well as the resistance tensor, |
967 |
|
|
\[ |
968 |
|
|
\left( {\begin{array}{*{20}c} |
969 |
|
|
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
970 |
|
|
0& 0.9270&-0.007063& 0.08585&0&0\\ |
971 |
|
|
0&-0.007063&0.7494&0.2057&0&0\\ |
972 |
|
|
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). |
973 |
|
|
\] |
974 |
|
|
where the units for translational, translation-rotation coupling and rotational |
975 |
|
|
tensors are $\frac{kcal \cdot fs}{mol \cdot \rm{\AA}^2}$, $\frac{kcal \cdot fs}{ |
976 |
|
|
mol \cdot \rm{\AA} \cdot rad}$ and $\frac{kcal \cdot fs}{mol \cdot rad^2}$ respe |
977 |
|
|
ctively. |
978 |
gezelter |
3299 |
|
979 |
gezelter |
3308 |
The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor) |
980 |
|
|
are essentially quantitative for translational diffusion of this model. |
981 |
|
|
Orientational correlation times under the Langevin rigid-body integrator |
982 |
|
|
are within 11\% of the values obtained from explicit solvent, but these |
983 |
|
|
models also exhibit some solvent inaccessible surface area in the |
984 |
|
|
explicitly-solvated case. |
985 |
|
|
|
986 |
gezelter |
3299 |
\subsubsection{Composite sphero-ellipsoids} |
987 |
|
|
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
988 |
gezelter |
3302 |
used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01} |
989 |
gezelter |
3299 |
|
990 |
gezelter |
3302 |
The diffusion constants and rotation relaxation times for |
991 |
xsun |
3298 |
different shaped molecules are shown in table \ref{tab:translation} |
992 |
|
|
and \ref{tab:rotation} to compare to the results calculated from NVE |
993 |
|
|
simulations. The theoretical values for sphere is calculated from the |
994 |
|
|
Stokes-Einstein law, the theoretical values for ellipsoid is |
995 |
gezelter |
3299 |
calculated from Perrin's fomula, The exact method is |
996 |
xsun |
3298 |
applied to the langevin dynamics simulations for sphere and ellipsoid, |
997 |
|
|
the bead model is applied to the simulation for dumbbell molecule, and |
998 |
|
|
the rough shell model is applied to ellipsoid, dumbbell, banana and |
999 |
|
|
lipid molecules. The results from all the langevin dynamics |
1000 |
|
|
simulations, including exact, bead model and rough shell, match the |
1001 |
|
|
theoretical values perfectly for all different shaped molecules. This |
1002 |
|
|
indicates that our simulation package for langevin dynamics is working |
1003 |
|
|
well. The approxiate methods ( bead model and rough shell model) are |
1004 |
|
|
accurate enough for the current simulations. The goal of the langevin |
1005 |
|
|
dynamics theory is to replace the explicit solvents by the friction |
1006 |
|
|
forces. We compared the dynamic properties of different shaped |
1007 |
|
|
molecules in langevin dynamics simulations with that in NVE |
1008 |
|
|
simulations. The results are reasonable close. Overall, the |
1009 |
|
|
translational diffusion constants calculated from langevin dynamics |
1010 |
|
|
simulations are very close to the values from the NVE simulation. For |
1011 |
|
|
sphere and lipid molecules, the diffusion constants are a little bit |
1012 |
|
|
off from the NVE simulation results. One possible reason is that the |
1013 |
|
|
calculation of the viscosity is very difficult to be accurate. Another |
1014 |
|
|
possible reason is that although we save very frequently during the |
1015 |
|
|
NVE simulations and run pretty long time simulations, there is only |
1016 |
|
|
one solute molecule in the system which makes the calculation for the |
1017 |
|
|
diffusion constant difficult. The sphere molecule behaves as a free |
1018 |
|
|
rotor in the solvent, so there is no rotation relaxation time |
1019 |
|
|
calculated from NVE simulations. The rotation relaxation time is not |
1020 |
|
|
very close to the NVE simulations results. The banana and lipid |
1021 |
|
|
molecules match the NVE simulations results pretty well. The mismatch |
1022 |
|
|
between langevin dynamics and NVE simulation for ellipsoid is possibly |
1023 |
|
|
caused by the slip boundary condition. For dumbbell, the mismatch is |
1024 |
|
|
caused by the size of the solvent molecule is pretty large compared to |
1025 |
|
|
dumbbell molecule in NVE simulations. |
1026 |
|
|
|
1027 |
|
|
According to our simulations, the langevin dynamics is a reliable |
1028 |
|
|
theory to apply to replace the explicit solvents, especially for the |
1029 |
|
|
translation properties. For large molecules, the rotation properties |
1030 |
|
|
are also mimiced reasonablly well. |
1031 |
|
|
|
1032 |
|
|
\begin{table*} |
1033 |
|
|
\begin{minipage}{\linewidth} |
1034 |
|
|
\begin{center} |
1035 |
gezelter |
3305 |
\caption{Translational diffusion constants (D) for the model systems |
1036 |
|
|
calculated using microcanonical simulations (with explicit solvent), |
1037 |
|
|
theoretical predictions, and Langevin simulations (with implicit solvent). |
1038 |
|
|
Analytical solutions for the exactly-solved hydrodynamics models are |
1039 |
|
|
from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936} |
1040 |
|
|
(ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq} |
1041 |
|
|
(dumbbell). The other model systems have no known analytic solution. |
1042 |
|
|
All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= |
1043 |
|
|
$10^{-4}$ \AA$^2$ / fs). } |
1044 |
|
|
\begin{tabular}{lccccccc} |
1045 |
xsun |
3298 |
\hline |
1046 |
gezelter |
3305 |
& \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\ |
1047 |
|
|
\cline{2-3} \cline{5-7} |
1048 |
|
|
model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\ |
1049 |
xsun |
3298 |
\hline |
1050 |
xsun |
3309 |
sphere & 0.261 & ? & & 2.59 & exact & 2.59 & 2.56 \\ |
1051 |
gezelter |
3305 |
ellipsoid & 0.255 & 2.44 & & 2.34 & exact & 2.34 & 2.37 \\ |
1052 |
|
|
& 0.255 & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\ |
1053 |
xsun |
3309 |
dumbbell & 0.322 & ? & & 1.57 & bead model & 1.57 & 1.57 \\ |
1054 |
|
|
& 0.322 & ? & & 1.57 & rough shell & ? & ? \\ |
1055 |
gezelter |
3305 |
banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\ |
1056 |
|
|
lipid & 0.349 & 0.96 & & & rough shell & 1.33 & 1.32 \\ |
1057 |
xsun |
3298 |
\end{tabular} |
1058 |
|
|
\label{tab:translation} |
1059 |
|
|
\end{center} |
1060 |
|
|
\end{minipage} |
1061 |
|
|
\end{table*} |
1062 |
|
|
|
1063 |
|
|
\begin{table*} |
1064 |
|
|
\begin{minipage}{\linewidth} |
1065 |
|
|
\begin{center} |
1066 |
gezelter |
3305 |
\caption{Orientational relaxation times ($\tau$) for the model systems using |
1067 |
|
|
microcanonical simulation (with explicit solvent), theoretical |
1068 |
|
|
predictions, and Langevin simulations (with implicit solvent). All |
1069 |
|
|
relaxation times are for the rotational correlation function with |
1070 |
|
|
$\ell = 2$ and are reported in units of ps. The ellipsoidal model has |
1071 |
|
|
an exact solution for the orientational correlation time due to |
1072 |
|
|
Perrin, but the other model systems have no known analytic solution.} |
1073 |
|
|
\begin{tabular}{lccccccc} |
1074 |
xsun |
3298 |
\hline |
1075 |
gezelter |
3305 |
& \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\ |
1076 |
|
|
\cline{2-3} \cline{5-7} |
1077 |
|
|
model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic & simulation \\ |
1078 |
xsun |
3298 |
\hline |
1079 |
xsun |
3309 |
sphere & 0.261 & & & 9.06 & exact & 9.06 & 9.11 \\ |
1080 |
gezelter |
3305 |
ellipsoid & 0.255 & 46.7 & & 22.0 & exact & 22.0 & 22.2 \\ |
1081 |
|
|
& 0.255 & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\ |
1082 |
xsun |
3309 |
dumbbell & 0.322 & 14.0 & & & bead model & 52.3 & 52.8 \\ |
1083 |
|
|
& 0.322 & 14.0 & & & rough shell & ? & ? \\ |
1084 |
gezelter |
3305 |
banana & 0.298 & 63.8 & & & rough shell & 70.9 & 70.9 \\ |
1085 |
|
|
lipid & 0.349 & 78.0 & & & rough shell & 76.9 & 77.9 \\ |
1086 |
|
|
\hline |
1087 |
xsun |
3298 |
\end{tabular} |
1088 |
|
|
\label{tab:rotation} |
1089 |
|
|
\end{center} |
1090 |
|
|
\end{minipage} |
1091 |
|
|
\end{table*} |
1092 |
|
|
|
1093 |
|
|
Langevin dynamics simulations are applied to study the formation of |
1094 |
|
|
the ripple phase of lipid membranes. The initial configuration is |
1095 |
|
|
taken from our molecular dynamics studies on lipid bilayers with |
1096 |
|
|
lennard-Jones sphere solvents. The solvent molecules are excluded from |
1097 |
|
|
the system, the experimental value of water viscosity is applied to |
1098 |
|
|
mimic the heat bath. Fig. XXX is the snapshot of the stable |
1099 |
|
|
configuration of the system, the ripple structure stayed stable after |
1100 |
|
|
100 ns run. The efficiency of the simulation is increased by one order |
1101 |
|
|
of magnitude. |
1102 |
|
|
|
1103 |
tim |
2746 |
\section{Conclusions} |
1104 |
|
|
|
1105 |
tim |
2999 |
We have presented a new Langevin algorithm by incorporating the |
1106 |
|
|
hydrodynamics properties of arbitrary shaped molecules into an |
1107 |
gezelter |
3308 |
advanced symplectic integration scheme. Further studies in systems |
1108 |
|
|
involving banana shaped molecules illustrated that the dynamic |
1109 |
|
|
properties could be preserved by using this new algorithm as an |
1110 |
|
|
implicit solvent model. |
1111 |
tim |
2999 |
|
1112 |
|
|
|
1113 |
tim |
2746 |
\section{Acknowledgments} |
1114 |
|
|
Support for this project was provided by the National Science |
1115 |
|
|
Foundation under grant CHE-0134881. T.L. also acknowledges the |
1116 |
|
|
financial support from center of applied mathematics at University |
1117 |
|
|
of Notre Dame. |
1118 |
|
|
\newpage |
1119 |
|
|
|
1120 |
gezelter |
3305 |
\bibliographystyle{jcp} |
1121 |
tim |
2746 |
\bibliography{langevin} |
1122 |
|
|
|
1123 |
|
|
\end{document} |