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