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