1 |
|
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
2 |
< |
%\documentclass[aps,prb,preprint]{revtex4} |
2 |
> |
%\documentclass[aps,prb,preprint,amsmath,amssymb,endfloats]{revtex4} |
3 |
|
\documentclass[11pt]{article} |
4 |
– |
\usepackage{endfloat} |
4 |
|
\usepackage{amsmath} |
5 |
|
\usepackage{amssymb} |
6 |
|
\usepackage{times} |
7 |
|
\usepackage{mathptmx} |
8 |
|
\usepackage{setspace} |
9 |
< |
\usepackage{tabularx} |
9 |
> |
\usepackage{endfloat} |
10 |
> |
%\usepackage{tabularx} |
11 |
|
\usepackage{graphicx} |
12 |
< |
\usepackage{booktabs} |
13 |
< |
\usepackage{bibentry} |
14 |
< |
\usepackage{mathrsfs} |
12 |
> |
%\usepackage{booktabs} |
13 |
> |
%\usepackage{bibentry} |
14 |
> |
%\usepackage{mathrsfs} |
15 |
|
\usepackage[ref]{overcite} |
16 |
|
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm |
17 |
|
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight |
18 |
|
9.0in \textwidth 6.5in \brokenpenalty=10000 |
19 |
< |
\renewcommand{\baselinestretch}{1.2} |
19 |
> |
|
20 |
|
\renewcommand\citemid{\ } % no comma in optional referenc note |
21 |
|
|
22 |
|
\begin{document} |
23 |
|
|
24 |
< |
\title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape } |
24 |
> |
\title{Langevin dynamics for rigid bodies of arbitrary shape} |
25 |
|
|
26 |
< |
\author{Xiuquan Sun, Teng Lin and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
27 |
< |
gezelter@nd.edu} \\ |
28 |
< |
Department of Chemistry and Biochemistry\\ |
26 |
> |
\author{Xiuquan Sun, Teng Lin and J. Daniel |
27 |
> |
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ |
28 |
> |
Department of Chemistry and Biochemistry,\\ |
29 |
|
University of Notre Dame\\ |
30 |
|
Notre Dame, Indiana 46556} |
31 |
|
|
32 |
|
\date{\today} |
33 |
|
|
34 |
– |
\maketitle \doublespacing |
34 |
|
|
35 |
< |
\begin{abstract} |
35 |
> |
\maketitle |
36 |
|
|
37 |
+ |
|
38 |
+ |
|
39 |
+ |
\begin{abstract} |
40 |
+ |
We present an algorithm for carrying out Langevin dynamics simulations |
41 |
+ |
on complex rigid bodies by incorporating the hydrodynamic resistance |
42 |
+ |
tensors for arbitrary shapes into an advanced symplectic integration |
43 |
+ |
scheme. The integrator gives quantitative agreement with both |
44 |
+ |
analytic and approximate hydrodynamic theories for a number of model |
45 |
+ |
rigid bodies, and works well at reproducing the solute dynamical |
46 |
+ |
properties (diffusion constants, and orientational relaxation times) |
47 |
+ |
obtained from explicitly-solvated simulations. |
48 |
|
\end{abstract} |
49 |
|
|
50 |
|
\newpage |
51 |
|
|
52 |
+ |
|
53 |
+ |
|
54 |
|
%\narrowtext |
55 |
|
|
56 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
57 |
|
% BODY OF TEXT |
58 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
59 |
|
|
60 |
+ |
\begin{doublespace} |
61 |
+ |
|
62 |
|
\section{Introduction} |
63 |
|
|
64 |
|
%applications of langevin dynamics |
65 |
< |
Langevin dynamics, which mimics a simple heat bath with stochastic and |
65 |
> |
Langevin dynamics, which mimics a heat bath using both stochastic and |
66 |
|
dissipative forces, has been applied in a variety of situations as an |
67 |
|
alternative to molecular dynamics with explicit solvent molecules. |
68 |
|
The stochastic treatment of the solvent allows the use of simulations |
69 |
< |
with substantially longer time and length scales. In general, the |
69 |
> |
with substantially longer time and length scales. In general, the |
70 |
|
dynamic and structural properties obtained from Langevin simulations |
71 |
|
agree quite well with similar properties obtained from explicit |
72 |
|
solvent simulations. |
75 |
|
study of met-enkephalin in which Langevin simulations predicted |
76 |
|
dynamical properties that were largely in agreement with explicit |
77 |
|
solvent simulations.\cite{Shen2002} By applying Langevin dynamics with |
78 |
< |
the UNRES model, Liow and his coworkers suggest that protein folding |
78 |
> |
the UNRES model, Liwo and his coworkers suggest that protein folding |
79 |
|
pathways can be explored within a reasonable amount of |
80 |
|
time.\cite{Liwo2005} |
81 |
|
|
90 |
|
folding/unfolding studies and discovered a higher free energy barrier |
91 |
|
between the native and denatured states.\cite{HuseyinKaya07012005} |
92 |
|
|
93 |
< |
Because of its stability against noise, Langevin dynamics has also |
80 |
< |
proven useful for studying remagnetization processes in various |
81 |
< |
systems.\cite{Palacios1998,Berkov2002,Denisov2003} [Check: For |
82 |
< |
instance, the oscillation power spectrum of nanoparticles from |
83 |
< |
Langevin dynamics has the same peak frequencies for different wave |
84 |
< |
vectors, which recovers the property of magnetic excitations in small |
85 |
< |
finite structures.\cite{Berkov2005a}] |
86 |
< |
|
87 |
< |
In typical LD simulations, the friction and random forces on |
93 |
> |
In typical LD simulations, the friction and random ($f_r$) forces on |
94 |
|
individual atoms are taken from Stokes' law, |
95 |
|
\begin{eqnarray} |
96 |
< |
m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\ |
97 |
< |
\langle R(t) \rangle & = & 0 \\ |
98 |
< |
\langle R(t) R(t') \rangle & = & 2 k_B T \xi m \delta(t - t') |
96 |
> |
m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + f_r(t) \notag \\ |
97 |
> |
\langle f_r(t) \rangle & = & 0 \\ |
98 |
> |
\langle f_r(t) f_r(t') \rangle & = & 2 k_B T \xi m \delta(t - t') \notag |
99 |
|
\end{eqnarray} |
100 |
< |
where $\xi \approx 6 \pi \eta a$. Here $\eta$ is the viscosity of the |
101 |
< |
implicit solvent, and $a$ is the hydrodynamic radius of the atom. |
100 |
> |
where $\xi \approx 6 \pi \eta \rho$. Here $\eta$ is the viscosity of the |
101 |
> |
implicit solvent, and $\rho$ is the hydrodynamic radius of the atom. |
102 |
|
|
103 |
|
The use of rigid substructures,\cite{Chun:2000fj} |
104 |
< |
coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunGezelter08} |
105 |
< |
and ellipsoidal representations of protein side chains~\cite{Fogolari:1996lr} |
106 |
< |
has made the use of the Stokes-Einstein approximation problematic. A |
107 |
< |
rigid substructure moves as a single unit with orientational as well |
108 |
< |
as translational degrees of freedom. This requires a more general |
109 |
< |
treatment of the hydrodynamics than the spherical approximation |
110 |
< |
provides. The atoms involved in a rigid or coarse-grained structure |
111 |
< |
should properly have solvent-mediated interactions with each |
112 |
< |
other. The theory of interactions {\it between} bodies moving through |
113 |
< |
a fluid has been developed over the past century and has been applied |
114 |
< |
to simulations of Brownian |
115 |
< |
motion.\cite{FIXMAN:1986lr,Ramachandran1996} |
104 |
> |
coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunX._jp0762020} |
105 |
> |
and ellipsoidal representations of protein side |
106 |
> |
chains~\cite{Fogolari:1996lr} has made the use of the Stokes-Einstein |
107 |
> |
approximation problematic. A rigid substructure moves as a single |
108 |
> |
unit with orientational as well as translational degrees of freedom. |
109 |
> |
This requires a more general treatment of the hydrodynamics than the |
110 |
> |
spherical approximation provides. Also, the atoms involved in a rigid |
111 |
> |
or coarse-grained structure have solvent-mediated interactions with |
112 |
> |
each other, and these interactions are ignored if all atoms are |
113 |
> |
treated as separate spherical particles. The theory of interactions |
114 |
> |
{\it between} bodies moving through a fluid has been developed over |
115 |
> |
the past century and has been applied to simulations of Brownian |
116 |
> |
motion.\cite{FIXMAN:1986lr,Ramachandran1996} |
117 |
|
|
118 |
< |
In order to account for the diffusion anisotropy of arbitrarily-shaped |
119 |
< |
particles, Fernandes and Garc\'{i}a de la Torre improved the original |
120 |
< |
Brownian dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by |
118 |
> |
In order to account for the diffusion anisotropy of complex shapes, |
119 |
> |
Fernandes and Garc\'{i}a de la Torre improved an earlier Brownian |
120 |
> |
dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by |
121 |
|
incorporating a generalized $6\times6$ diffusion tensor and |
122 |
|
introducing a rotational evolution scheme consisting of three |
123 |
|
consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are |
128 |
|
assumption of zero average acceleration is not always true for |
129 |
|
cooperative motion which is common in proteins. An inertial Brownian |
130 |
|
dynamics (IBD) was proposed to address this issue by adding an |
131 |
< |
inertial correction term.\cite{Beard2000} As a complement to IBD which |
132 |
< |
has a lower bound in time step because of the inertial relaxation |
133 |
< |
time, long-time-step inertial dynamics (LTID) can be used to |
134 |
< |
investigate the inertial behavior of linked polymer segments in a low |
135 |
< |
friction regime.\cite{Beard2000} LTID can also deal with the |
131 |
> |
inertial correction term.\cite{Beard2000} As a complement to IBD, |
132 |
> |
which has a lower bound in time step because of the inertial |
133 |
> |
relaxation time, long-time-step inertial dynamics (LTID) can be used |
134 |
> |
to investigate the inertial behavior of linked polymer segments in a |
135 |
> |
low friction regime.\cite{Beard2000} LTID can also deal with the |
136 |
|
rotational dynamics for nonskew bodies without translation-rotation |
137 |
|
coupling by separating the translation and rotation motion and taking |
138 |
< |
advantage of the analytical solution of hydrodynamics |
138 |
> |
advantage of the analytical solution of hydrodynamic |
139 |
|
properties. However, typical nonskew bodies like cylinders and |
140 |
|
ellipsoids are inadequate to represent most complex macromolecular |
141 |
< |
assemblies. There is therefore a need for incorporating the |
142 |
< |
hydrodynamics of complex (and potentially skew) rigid bodies in the |
143 |
< |
library of methods available for performing Langevin simulations. |
141 |
> |
assemblies. Therefore, the goal of this work is to adapt some of the |
142 |
> |
hydrodynamic methodologies developed to treat Brownian motion of |
143 |
> |
complex assemblies into a Langevin integrator for rigid bodies with |
144 |
> |
arbitrary shapes. |
145 |
|
|
146 |
|
\subsection{Rigid Body Dynamics} |
147 |
|
Rigid bodies are frequently involved in the modeling of large |
148 |
|
collections of particles that move as a single unit. In molecular |
149 |
|
simulations, rigid bodies have been used to simplify protein-protein |
150 |
|
docking,\cite{Gray2003} and lipid bilayer |
151 |
< |
simulations.\cite{SunGezelter08} Many of the water models in common |
151 |
> |
simulations.\cite{SunX._jp0762020} Many of the water models in common |
152 |
|
use are also rigid-body |
153 |
|
models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are |
154 |
< |
typically evolved using constraints rather than rigid body equations |
155 |
< |
of motion. |
154 |
> |
typically evolved in molecular dynamics simulations using constraints |
155 |
> |
rather than rigid body equations of motion. |
156 |
|
|
157 |
|
Euler angles are a natural choice to describe the rotational degrees |
158 |
|
of freedom. However, due to $\frac{1}{\sin \theta}$ singularities, the |
177 |
|
In order to develop a stable and efficient integration scheme that |
178 |
|
preserves most constants of the motion, symplectic propagators are |
179 |
|
necessary. By introducing a conjugate momentum to the rotation matrix |
180 |
< |
$Q$ and re-formulating Hamilton's equations, a symplectic |
180 |
> |
${\bf Q}$ and re-formulating Hamilton's equations, a symplectic |
181 |
|
orientational integrator, RSHAKE,\cite{Kol1997} was proposed to evolve |
182 |
|
rigid bodies on a constraint manifold by iteratively satisfying the |
183 |
< |
orthogonality constraint $Q^T Q = 1$. An alternative method using the |
184 |
< |
quaternion representation was developed by Omelyan.\cite{Omelyan1998} |
185 |
< |
However, both of these methods are iterative and suffer from some |
186 |
< |
related inefficiencies. A symplectic Lie-Poisson integrator for rigid |
187 |
< |
bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes |
188 |
< |
most of the limitations mentioned above and is therefore the basis for |
189 |
< |
our Langevin integrator. |
183 |
> |
orthogonality constraint ${\bf Q}^T {\bf Q} = 1$. An alternative |
184 |
> |
method using the quaternion representation was developed by |
185 |
> |
Omelyan.\cite{Omelyan1998} However, both of these methods are |
186 |
> |
iterative and suffer from some related inefficiencies. A symplectic |
187 |
> |
Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et |
188 |
> |
al.}\cite{Dullweber1997} removes most of the limitations mentioned |
189 |
> |
above and is therefore the basis for our Langevin integrator. |
190 |
|
|
191 |
|
The goal of the present work is to develop a Langevin dynamics |
192 |
< |
algorithm for arbitrary-shaped rigid particles by integrating the |
193 |
< |
accurate estimation of friction tensor from hydrodynamics theory into |
194 |
< |
a symplectic rigid body dynamics propagator. In the sections below, |
195 |
< |
we review some of the theory of hydrodynamic tensors developed |
192 |
> |
algorithm for arbitrary-shaped rigid particles by integrating an |
193 |
> |
accurate estimate of the friction tensor from hydrodynamics theory |
194 |
> |
into a symplectic rigid body dynamics propagator. In the sections |
195 |
> |
below, we review some of the theory of hydrodynamic tensors developed |
196 |
|
primarily for Brownian simulations of multi-particle systems, we then |
197 |
|
present our integration method for a set of generalized Langevin |
198 |
|
equations of motion, and we compare the behavior of the new Langevin |
200 |
|
molecular dynamics. |
201 |
|
|
202 |
|
\subsection{\label{introSection:frictionTensor}The Friction Tensor} |
203 |
< |
Theoretically, a complete friction kernel can be determined using the |
204 |
< |
velocity autocorrelation function. However, this approach becomes |
205 |
< |
impractical when the solute becomes complex. Instead, various |
206 |
< |
approaches based on hydrodynamics have been developed to calculate the |
207 |
< |
friction coefficients. In general, the friction tensor $\Xi$ is a |
208 |
< |
$6\times 6$ matrix given by |
203 |
> |
Theoretically, a complete friction kernel for a solute particle can be |
204 |
> |
determined using the velocity autocorrelation function from a |
205 |
> |
simulation with explicit solvent molecules. However, this approach |
206 |
> |
becomes impractical when the solute becomes complex. Instead, various |
207 |
> |
approaches based on hydrodynamics have been developed to calculate |
208 |
> |
static friction coefficients. In general, the friction tensor $\Xi$ is |
209 |
> |
a $6\times 6$ matrix given by |
210 |
|
\begin{equation} |
211 |
|
\Xi = \left( \begin{array}{*{20}c} |
212 |
|
\Xi^{tt} & \Xi^{rt} \\ |
217 |
|
rotational resistance (friction) tensors respectively, while |
218 |
|
$\Xi^{tr}$ is translation-rotation coupling tensor and $\Xi^{rt}$ is |
219 |
|
rotation-translation coupling tensor. When a particle moves in a |
220 |
< |
fluid, it may experience friction force ($\mathbf{f}_f$) and torque |
221 |
< |
($\mathbf{\tau}_f$) in opposition to the directions of the velocity |
222 |
< |
($\mathbf{v}$) and body-fixed angular velocity ($\mathbf{\omega}$), |
220 |
> |
fluid, it may experience a friction force ($\mathbf{f}_f$) and torque |
221 |
> |
($\mathbf{\tau}_f$) in opposition to the velocity ($\mathbf{v}$) and |
222 |
> |
body-fixed angular velocity ($\mathbf{\omega}$), |
223 |
|
\begin{equation} |
224 |
|
\left( \begin{array}{l} |
225 |
|
\mathbf{f}_f \\ |
234 |
|
\end{equation} |
235 |
|
|
236 |
|
\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} |
237 |
< |
For a spherical particle under ``stick'' boundary conditions, the |
238 |
< |
translational and rotational friction tensors can be calculated from |
239 |
< |
Stokes' law, |
237 |
> |
For a spherical body under ``stick'' boundary conditions, |
238 |
> |
the translational and rotational friction tensors can be calculated |
239 |
> |
from Stokes' law, |
240 |
|
\begin{equation} |
241 |
+ |
\label{eq:StokesTranslation} |
242 |
|
\Xi^{tt} = \left( \begin{array}{*{20}c} |
243 |
< |
{6\pi \eta R} & 0 & 0 \\ |
244 |
< |
0 & {6\pi \eta R} & 0 \\ |
245 |
< |
0 & 0 & {6\pi \eta R} \\ |
243 |
> |
{6\pi \eta \rho} & 0 & 0 \\ |
244 |
> |
0 & {6\pi \eta \rho} & 0 \\ |
245 |
> |
0 & 0 & {6\pi \eta \rho} \\ |
246 |
|
\end{array} \right) |
247 |
|
\end{equation} |
248 |
|
and |
249 |
|
\begin{equation} |
250 |
+ |
\label{eq:StokesRotation} |
251 |
|
\Xi^{rr} = \left( \begin{array}{*{20}c} |
252 |
< |
{8\pi \eta R^3 } & 0 & 0 \\ |
253 |
< |
0 & {8\pi \eta R^3 } & 0 \\ |
254 |
< |
0 & 0 & {8\pi \eta R^3 } \\ |
252 |
> |
{8\pi \eta \rho^3 } & 0 & 0 \\ |
253 |
> |
0 & {8\pi \eta \rho^3 } & 0 \\ |
254 |
> |
0 & 0 & {8\pi \eta \rho^3 } \\ |
255 |
|
\end{array} \right) |
256 |
|
\end{equation} |
257 |
< |
where $\eta$ is the viscosity of the solvent and $R$ is the |
258 |
< |
hydrodynamic radius. |
257 |
> |
where $\eta$ is the viscosity of the solvent and $\rho$ is the |
258 |
> |
hydrodynamic radius. The presence of the rotational resistance tensor |
259 |
> |
implies that the spherical body has internal structure and |
260 |
> |
orientational degrees of freedom that must be propagated in time. For |
261 |
> |
non-structured spherical bodies (i.e. the atoms in a traditional |
262 |
> |
molecular dynamics simulation) these degrees of freedom do not exist. |
263 |
|
|
264 |
|
Other non-spherical shapes, such as cylinders and ellipsoids, are |
265 |
< |
widely used as references for developing new hydrodynamics theories, |
265 |
> |
widely used as references for developing new hydrodynamic theories, |
266 |
|
because their properties can be calculated exactly. In 1936, Perrin |
267 |
< |
extended Stokes' law to general ellipsoids which are given in |
268 |
< |
Cartesian coordinates by~\cite{Perrin1934,Perrin1936} |
267 |
> |
extended Stokes' law to general |
268 |
> |
ellipsoids,\cite{Perrin1934,Perrin1936} described in Cartesian |
269 |
> |
coordinates as |
270 |
|
\begin{equation} |
271 |
|
\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1. |
272 |
|
\end{equation} |
273 |
|
Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the |
274 |
|
complexity of the elliptic integral, only uniaxial ellipsoids, either |
275 |
< |
prolate ($a \ge b = c$) or oblate ($a < b = c$), can be solved |
275 |
> |
prolate ($a \ge b = c$) or oblate ($a < b = c$), were solved |
276 |
|
exactly. Introducing an elliptic integral parameter $S$ for prolate, |
277 |
|
\begin{equation} |
278 |
|
S = \frac{2}{\sqrt{a^2 - b^2}} \ln \frac{a + \sqrt{a^2 - b^2}}{b}, |
281 |
|
\begin{equation} |
282 |
|
S = \frac{2}{\sqrt {b^2 - a^2 }} \arctan \frac{\sqrt {b^2 - a^2}}{a}, |
283 |
|
\end{equation} |
284 |
< |
ellipsoids, one can write down the translational and rotational |
285 |
< |
resistance tensors: |
286 |
< |
\begin{eqnarray*} |
284 |
> |
ellipsoids, it is possible to write down exact solutions for the |
285 |
> |
resistance tensors. As is the case for spherical bodies, the translational, |
286 |
> |
\begin{eqnarray} |
287 |
|
\Xi_a^{tt} & = & 16\pi \eta \frac{a^2 - b^2}{(2a^2 - b^2 )S - 2a}. \\ |
288 |
|
\Xi_b^{tt} = \Xi_c^{tt} & = & 32\pi \eta \frac{a^2 - b^2 }{(2a^2 - 3b^2 )S + 2a}, |
289 |
< |
\end{eqnarray*} |
290 |
< |
for oblate, and |
291 |
< |
\begin{eqnarray*} |
289 |
> |
\end{eqnarray} |
290 |
> |
and rotational, |
291 |
> |
\begin{eqnarray} |
292 |
|
\Xi_a^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^2 - b^2 )b^2}{2a - b^2 S}, \\ |
293 |
|
\Xi_b^{rr} = \Xi_c^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^4 - b^4)}{(2a^2 - b^2 )S - 2a} |
294 |
< |
\end{eqnarray*} |
295 |
< |
for prolate ellipsoids. For both spherical and ellipsoidal particles, |
296 |
< |
the translation-rotation and rotation-translation coupling tensors are |
297 |
< |
zero. |
294 |
> |
\end{eqnarray} |
295 |
> |
resistance tensors are diagonal $3 \times 3$ matrices. For both |
296 |
> |
spherical and ellipsoidal particles, the translation-rotation and |
297 |
> |
rotation-translation coupling tensors are zero. |
298 |
|
|
299 |
|
\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}} |
300 |
< |
Unlike spherical and other simply shaped molecules, there is no |
301 |
< |
analytical solution for the friction tensor for arbitrarily shaped |
302 |
< |
rigid molecules. The ellipsoid of revolution model and general |
287 |
< |
triaxial ellipsoid model have been used to approximate the |
300 |
> |
There is no analytical solution for the friction tensor for rigid |
301 |
> |
molecules of arbitrary shape. The ellipsoid of revolution and general |
302 |
> |
triaxial ellipsoid models have been widely used to approximate the |
303 |
|
hydrodynamic properties of rigid bodies. However, the mapping from all |
304 |
< |
possible ellipsoidal spaces, $r$-space, to all possible combination of |
305 |
< |
rotational diffusion coefficients, $D$-space, is not |
304 |
> |
possible ellipsoidal spaces ($r$-space) to all possible combinations |
305 |
> |
of rotational diffusion coefficients ($D$-space) is not |
306 |
|
unique.\cite{Wegener1979} Additionally, because there is intrinsic |
307 |
< |
coupling between translational and rotational motion of rigid bodies, |
308 |
< |
general ellipsoids are not always suitable for modeling arbitrarily |
309 |
< |
shaped rigid molecules. A number of studies have been devoted to |
310 |
< |
determining the friction tensor for irregularly shaped rigid bodies |
311 |
< |
using more advanced methods where the molecule of interest was modeled |
312 |
< |
by a combinations of spheres\cite{Carrasco1999} and the hydrodynamics |
313 |
< |
properties of the molecule can be calculated using the hydrodynamic |
314 |
< |
interaction tensor. |
307 |
> |
coupling between translational and rotational motion of {\it skew} |
308 |
> |
rigid bodies, general ellipsoids are not always suitable for modeling |
309 |
> |
rigid molecules. A number of studies have been devoted to determining |
310 |
> |
the friction tensor for irregular shapes using methods in which the |
311 |
> |
molecule of interest is modeled with a combination of |
312 |
> |
spheres\cite{Carrasco1999} and the hydrodynamic properties of the |
313 |
> |
molecule are then calculated using a set of two-point interaction |
314 |
> |
tensors. We have found the {\it bead} and {\it rough shell} models of |
315 |
> |
Carrasco and Garc\'{i}a de la Torre to be the most useful of these |
316 |
> |
methods,\cite{Carrasco1999} and we review the basic outline of the |
317 |
> |
rough shell approach here. A more thorough explanation can be found |
318 |
> |
in Ref. \citen{Carrasco1999}. |
319 |
|
|
320 |
< |
Consider a rigid assembly of $N$ beads immersed in a continuous |
321 |
< |
medium. Due to hydrodynamic interaction, the ``net'' velocity of $i$th |
322 |
< |
bead, $v'_i$ is different than its unperturbed velocity $v_i$, |
320 |
> |
Consider a rigid assembly of $N$ small beads moving through a |
321 |
> |
continuous medium. Due to hydrodynamic interactions between the |
322 |
> |
beads, the net velocity of the $i^\mathrm{th}$ bead relative to the |
323 |
> |
medium, ${\bf v}'_i$, is different than its unperturbed velocity ${\bf |
324 |
> |
v}_i$, |
325 |
|
\begin{equation} |
326 |
< |
v'_i = v_i - \sum\limits_{j \ne i} {T_{ij} F_j } |
326 |
> |
{\bf v}'_i = {\bf v}_i - \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j } |
327 |
|
\end{equation} |
328 |
< |
where $F_i$ is the frictional force, and $T_{ij}$ is the hydrodynamic |
329 |
< |
interaction tensor. The frictional force on the $i^\mathrm{th}$ bead |
330 |
< |
is proportional to its ``net'' velocity |
328 |
> |
where ${\bf F}_j$ is the frictional force on the medium due to bead $j$, and |
329 |
> |
${\bf T}_{ij}$ is the hydrodynamic interaction tensor between the two beads. |
330 |
> |
The frictional force felt by the $i^\mathrm{th}$ bead is proportional to |
331 |
> |
its net velocity |
332 |
|
\begin{equation} |
333 |
< |
F_i = \zeta _i v_i - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }. |
333 |
> |
{\bf F}_i = \xi_i {\bf v}_i - \xi_i \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }. |
334 |
|
\label{introEquation:tensorExpression} |
335 |
|
\end{equation} |
336 |
< |
This equation is the basis for deriving the hydrodynamic tensor. In |
337 |
< |
1930, Oseen and Burgers gave a simple solution to |
338 |
< |
Eq.~\ref{introEquation:tensorExpression} |
339 |
< |
\begin{equation} |
340 |
< |
T_{ij} = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij} |
341 |
< |
R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor} |
320 |
< |
\end{equation} |
321 |
< |
Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$. |
322 |
< |
A second order expression for element of different size was |
323 |
< |
introduced by Rotne and Prager\cite{Rotne1969} and improved by |
336 |
> |
Eq. (\ref{introEquation:tensorExpression}) defines the two-point |
337 |
> |
hydrodynamic tensor, ${\bf T}_{ij}$. There have been many proposed |
338 |
> |
solutions to this equation, including the simple solution given by |
339 |
> |
Oseen and Burgers in 1930 for two beads of identical radius. A second |
340 |
> |
order expression for beads of different hydrodynamic radii was |
341 |
> |
introduced by Rotne and Prager,\cite{Rotne1969} and improved by |
342 |
|
Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977} |
343 |
|
\begin{equation} |
344 |
< |
T_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I + |
345 |
< |
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma |
346 |
< |
_i^2 + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} - |
347 |
< |
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. |
344 |
> |
{\bf T}_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {{\bf I} + |
345 |
> |
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\rho |
346 |
> |
_i^2 + \rho_j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} - |
347 |
> |
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. |
348 |
|
\label{introEquation:RPTensorNonOverlapped} |
349 |
|
\end{equation} |
350 |
< |
Both of the Eq.~\ref{introEquation:oseenTensor} and |
351 |
< |
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption |
352 |
< |
$R_{ij} \ge \sigma _i + \sigma _j$. An alternative expression for |
353 |
< |
overlapping beads with the same radius, $\sigma$, is given by |
350 |
> |
Here ${\bf R}_{ij}$ is the distance vector between beads $i$ and $j$. Both |
351 |
> |
the Oseen-Burgers tensor and |
352 |
> |
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption that |
353 |
> |
the beads do not overlap ($R_{ij} \ge \rho_i + \rho_j$). |
354 |
> |
|
355 |
> |
To calculate the resistance tensor for a body represented as the union |
356 |
> |
of many non-overlapping beads, we first pick an arbitrary origin $O$ |
357 |
> |
and then construct a $3N \times 3N$ supermatrix consisting of $N |
358 |
> |
\times N$ ${\bf B}_{ij}$ blocks |
359 |
|
\begin{equation} |
360 |
< |
T_{ij} = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 - |
361 |
< |
\frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I + |
362 |
< |
\frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right] |
363 |
< |
\label{introEquation:RPTensorOverlapped} |
360 |
> |
{\bf B} = \left( \begin{array}{*{20}c} |
361 |
> |
{\bf B}_{11} & \ldots & {\bf B}_{1N} \\ |
362 |
> |
\vdots & \ddots & \vdots \\ |
363 |
> |
{\bf B}_{N1} & \cdots & {\bf B}_{NN} |
364 |
> |
\end{array} \right) |
365 |
|
\end{equation} |
366 |
< |
To calculate the resistance tensor at an arbitrary origin $O$, we |
367 |
< |
construct a $3N \times 3N$ matrix consisting of $N \times N$ |
344 |
< |
$B_{ij}$ blocks |
366 |
> |
${\bf B}_{ij}$ is a version of the hydrodynamic tensor which includes the |
367 |
> |
self-contributions for spheres, |
368 |
|
\begin{equation} |
369 |
< |
B = \left( \begin{array}{*{20}c} |
370 |
< |
B_{11} & \ldots & B_{1N} \\ |
371 |
< |
\vdots & \ddots & \vdots \\ |
372 |
< |
B_{N1} & \cdots & B_{NN} \\ |
369 |
> |
{\bf B}_{ij} = \delta _{ij} \frac{{\bf I}}{{6\pi \eta R_{ij}}} + (1 - \delta_{ij} |
370 |
> |
){\bf T}_{ij} |
371 |
> |
\end{equation} |
372 |
> |
where $\delta_{ij}$ is the Kronecker delta function. Inverting the |
373 |
> |
${\bf B}$ matrix, we obtain |
374 |
> |
\begin{equation} |
375 |
> |
{\bf C} = {\bf B}^{ - 1} = \left(\begin{array}{*{20}c} |
376 |
> |
{\bf C}_{11} & \ldots & {\bf C}_{1N} \\ |
377 |
> |
\vdots & \ddots & \vdots \\ |
378 |
> |
{\bf C}_{N1} & \cdots & {\bf C}_{NN} |
379 |
|
\end{array} \right), |
380 |
|
\end{equation} |
381 |
< |
where $B_{ij}$ is given by |
381 |
> |
which can be partitioned into $N \times N$ blocks labeled ${\bf C}_{ij}$. |
382 |
> |
(Each of the ${\bf C}_{ij}$ blocks is a $3 \times 3$ matrix.) Using the |
383 |
> |
skew matrix, |
384 |
|
\begin{equation} |
385 |
< |
B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij} |
386 |
< |
)T_{ij} |
385 |
> |
{\bf U}_i = \left(\begin{array}{*{20}c} |
386 |
> |
0 & -z_i & y_i \\ |
387 |
> |
z_i & 0 & - x_i \\ |
388 |
> |
-y_i & x_i & 0 |
389 |
> |
\end{array}\right) |
390 |
> |
\label{eq:skewMatrix} |
391 |
|
\end{equation} |
357 |
– |
where $\delta _{ij}$ is the Kronecker delta function. Inverting the |
358 |
– |
$B$ matrix, we obtain |
359 |
– |
\[ |
360 |
– |
C = B^{ - 1} = \left( {\begin{array}{*{20}c} |
361 |
– |
{C_{11} } & \ldots & {C_{1N} } \\ |
362 |
– |
\vdots & \ddots & \vdots \\ |
363 |
– |
{C_{N1} } & \cdots & {C_{NN} } \\ |
364 |
– |
\end{array}} \right), |
365 |
– |
\] |
366 |
– |
which can be partitioned into $N \times N$ $3 \times 3$ block |
367 |
– |
$C_{ij}$. With the help of $C_{ij}$ and the skew matrix $U_i$ |
368 |
– |
\[ |
369 |
– |
U_i = \left( {\begin{array}{*{20}c} |
370 |
– |
0 & { - z_i } & {y_i } \\ |
371 |
– |
{z_i } & 0 & { - x_i } \\ |
372 |
– |
{ - y_i } & {x_i } & 0 \\ |
373 |
– |
\end{array}} \right) |
374 |
– |
\] |
392 |
|
where $x_i$, $y_i$, $z_i$ are the components of the vector joining |
393 |
< |
bead $i$ and origin $O$, the elements of resistance tensor at |
394 |
< |
arbitrary origin $O$ can be written as |
393 |
> |
bead $i$ and origin $O$, the elements of the resistance tensor (at the |
394 |
> |
arbitrary origin $O$) can be written as |
395 |
|
\begin{eqnarray} |
396 |
|
\label{introEquation:ResistanceTensorArbitraryOrigin} |
397 |
< |
\Xi _{}^{tt} & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\ |
398 |
< |
\Xi _{}^{tr} & = & \Xi _{}^{rt} = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\ |
399 |
< |
\Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } |
383 |
< |
U_j + 6 \eta V {\bf I}. \notag |
397 |
> |
\Xi^{tt} & = & \sum\limits_i {\sum\limits_j {{\bf C}_{ij} } } \notag , \\ |
398 |
> |
\Xi^{tr} = \Xi _{}^{rt} & = & \sum\limits_i {\sum\limits_j {{\bf U}_i {\bf C}_{ij} } } , \\ |
399 |
> |
\Xi^{rr} & = & -\sum\limits_i \sum\limits_j {\bf U}_i {\bf C}_{ij} {\bf U}_j + 6 \eta V {\bf I}. \notag |
400 |
|
\end{eqnarray} |
401 |
< |
The final term in the expression for $\Xi^{rr}$ is correction that |
402 |
< |
accounts for errors in the rotational motion of certain kinds of bead |
403 |
< |
models. The additive correction uses the solvent viscosity ($\eta$) |
404 |
< |
as well as the total volume of the beads that contribute to the |
389 |
< |
hydrodynamic model, |
401 |
> |
The final term in the expression for $\Xi^{rr}$ is a correction that |
402 |
> |
accounts for errors in the rotational motion of the bead models. The |
403 |
> |
additive correction uses the solvent viscosity ($\eta$) as well as the |
404 |
> |
total volume of the beads that contribute to the hydrodynamic model, |
405 |
|
\begin{equation} |
406 |
< |
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3, |
406 |
> |
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3, |
407 |
|
\end{equation} |
408 |
< |
where $\sigma_i$ is the radius of bead $i$. This correction term was |
408 |
> |
where $\rho_i$ is the radius of bead $i$. This correction term was |
409 |
|
rigorously tested and compared with the analytical results for |
410 |
< |
two-sphere and ellipsoidal systems by Garcia de la Torre and |
410 |
> |
two-sphere and ellipsoidal systems by Garc\'{i}a de la Torre and |
411 |
|
Rodes.\cite{Torre:1983lr} |
412 |
|
|
413 |
< |
|
414 |
< |
The resistance tensor depends on the origin to which they refer. The |
415 |
< |
proper location for applying the friction force is the center of |
416 |
< |
resistance (or center of reaction), at which the trace of rotational |
417 |
< |
resistance tensor, $ \Xi ^{rr}$ reaches a minimum value. |
418 |
< |
Mathematically, the center of resistance is defined as an unique |
419 |
< |
point of the rigid body at which the translation-rotation coupling |
405 |
< |
tensors are symmetric, |
413 |
> |
In general, resistance tensors depend on the origin at which they were |
414 |
> |
computed. However, the proper location for applying the friction |
415 |
> |
force is the center of resistance, the special point at which the |
416 |
> |
trace of rotational resistance tensor, $\Xi^{rr}$ reaches a minimum |
417 |
> |
value. Mathematically, the center of resistance can also be defined |
418 |
> |
as the unique point for a rigid body at which the translation-rotation |
419 |
> |
coupling tensors are symmetric, |
420 |
|
\begin{equation} |
421 |
< |
\Xi^{tr} = \left( {\Xi^{tr} } \right)^T |
421 |
> |
\Xi^{tr} = \left(\Xi^{tr}\right)^T |
422 |
|
\label{introEquation:definitionCR} |
423 |
|
\end{equation} |
424 |
< |
From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, |
425 |
< |
we can easily derive that the translational resistance tensor is |
426 |
< |
origin independent, while the rotational resistance tensor and |
424 |
> |
From Eq. \ref{introEquation:ResistanceTensorArbitraryOrigin}, we can |
425 |
> |
easily derive that the {\it translational} resistance tensor is origin |
426 |
> |
independent, while the rotational resistance tensor and |
427 |
|
translation-rotation coupling resistance tensor depend on the |
428 |
< |
origin. Given the resistance tensor at an arbitrary origin $O$, and |
429 |
< |
a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can |
430 |
< |
obtain the resistance tensor at $P$ by |
431 |
< |
\begin{equation} |
432 |
< |
\begin{array}{l} |
433 |
< |
\Xi _P^{tt} = \Xi _O^{tt} \\ |
434 |
< |
\Xi _P^{tr} = \Xi _P^{rt} = \Xi _O^{tr} - U_{OP} \Xi _O^{tt} \\ |
435 |
< |
\Xi _P^{rr} = \Xi _O^{rr} - U_{OP} \Xi _O^{tt} U_{OP} + \Xi _O^{tr} U_{OP} - U_{OP} \Xi _O^{{tr} ^{^T }} \\ |
436 |
< |
\end{array} |
437 |
< |
\label{introEquation:resistanceTensorTransformation} |
438 |
< |
\end{equation} |
439 |
< |
where |
440 |
< |
\[ |
441 |
< |
U_{OP} = \left( {\begin{array}{*{20}c} |
442 |
< |
0 & { - z_{OP} } & {y_{OP} } \\ |
443 |
< |
{z_{OP} } & 0 & { - x_{OP} } \\ |
444 |
< |
{ - y_{OP} } & {x_{OP} } & 0 \\ |
445 |
< |
\end{array}} \right) |
446 |
< |
\] |
447 |
< |
Using Eq.~\ref{introEquation:definitionCR} and |
448 |
< |
Eq.~\ref{introEquation:resistanceTensorTransformation}, one can |
449 |
< |
locate the position of center of resistance, |
450 |
< |
\begin{eqnarray*} |
451 |
< |
\left( \begin{array}{l} |
452 |
< |
x_{OR} \\ |
453 |
< |
y_{OR} \\ |
454 |
< |
z_{OR} \\ |
455 |
< |
\end{array} \right) & = &\left( \begin{array}{*{20}c} |
456 |
< |
{(\Xi _O^{rr} )_{yy} + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} } \\ |
457 |
< |
{ - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz} + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} } \\ |
458 |
< |
{ - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx} + (\Xi _O^{rr} )_{yy} } \\ |
459 |
< |
\end{array} \right)^{ - 1} \\ |
446 |
< |
& & \left( \begin{array}{l} |
447 |
< |
(\Xi _O^{tr} )_{yz} - (\Xi _O^{tr} )_{zy} \\ |
448 |
< |
(\Xi _O^{tr} )_{zx} - (\Xi _O^{tr} )_{xz} \\ |
449 |
< |
(\Xi _O^{tr} )_{xy} - (\Xi _O^{tr} )_{yx} \\ |
450 |
< |
\end{array} \right) \\ |
451 |
< |
\end{eqnarray*} |
428 |
> |
origin. Given the resistance tensor at an arbitrary origin $O$, and a |
429 |
> |
vector ,${\bf r}_{OP} = (x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we |
430 |
> |
can obtain the resistance tensor at $P$ by |
431 |
> |
\begin{eqnarray} |
432 |
> |
\label{introEquation:resistanceTensorTransformation} |
433 |
> |
\Xi_P^{tt} & = & \Xi_O^{tt} \notag \\ |
434 |
> |
\Xi_P^{tr} = \Xi_P^{rt} & = & \Xi_O^{tr} - {\bf U}_{OP} \Xi _O^{tt} \\ |
435 |
> |
\Xi_P^{rr} & = &\Xi_O^{rr} - {\bf U}_{OP} \Xi_O^{tt} {\bf U}_{OP} |
436 |
> |
+ \Xi_O^{tr} {\bf U}_{OP} - {\bf U}_{OP} \left( \Xi_O^{tr} |
437 |
> |
\right)^{^T} \notag |
438 |
> |
\end{eqnarray} |
439 |
> |
where ${\bf U}_{OP}$ is the skew matrix (Eq. (\ref{eq:skewMatrix})) |
440 |
> |
for the vector between the origin $O$ and the point $P$. Using |
441 |
> |
Eqs.~\ref{introEquation:definitionCR}~and~\ref{introEquation:resistanceTensorTransformation}, |
442 |
> |
one can locate the position of center of resistance, |
443 |
> |
\begin{equation*} |
444 |
> |
\left(\begin{array}{l} |
445 |
> |
x_{OR} \\ |
446 |
> |
y_{OR} \\ |
447 |
> |
z_{OR} |
448 |
> |
\end{array}\right) = |
449 |
> |
\left(\begin{array}{*{20}c} |
450 |
> |
(\Xi_O^{rr})_{yy} + (\Xi_O^{rr})_{zz} & -(\Xi_O^{rr})_{xy} & -(\Xi_O^{rr})_{xz} \\ |
451 |
> |
-(\Xi_O^{rr})_{xy} & (\Xi_O^{rr})_{zz} + (\Xi_O^{rr})_{xx} & -(\Xi_O^{rr})_{yz} \\ |
452 |
> |
-(\Xi_O^{rr})_{xz} & -(\Xi_O^{rr})_{yz} & (\Xi_O^{rr})_{xx} + (\Xi_O^{rr})_{yy} \\ |
453 |
> |
\end{array}\right)^{-1} |
454 |
> |
\left(\begin{array}{l} |
455 |
> |
(\Xi_O^{tr})_{yz} - (\Xi_O^{tr})_{zy} \\ |
456 |
> |
(\Xi_O^{tr})_{zx} - (\Xi_O^{tr})_{xz} \\ |
457 |
> |
(\Xi_O^{tr})_{xy} - (\Xi_O^{tr})_{yx} |
458 |
> |
\end{array}\right) |
459 |
> |
\end{equation*} |
460 |
|
where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector |
461 |
|
joining center of resistance $R$ and origin $O$. |
462 |
|
|
463 |
+ |
For a general rigid molecular substructure, finding the $6 \times 6$ |
464 |
+ |
resistance tensor can be a computationally demanding task. First, a |
465 |
+ |
lattice of small beads that extends well beyond the boundaries of the |
466 |
+ |
rigid substructure is created. The lattice is typically composed of |
467 |
+ |
0.25 \AA\ beads on a dense FCC lattice. The lattice constant is taken |
468 |
+ |
to be the bead diameter, so that adjacent beads are touching, but do |
469 |
+ |
not overlap. To make a shape corresponding to the rigid structure, |
470 |
+ |
beads that sit on lattice sites that are outside the van der Waals |
471 |
+ |
radii of all of the atoms comprising the rigid body are excluded from |
472 |
+ |
the calculation. |
473 |
|
|
474 |
+ |
For large structures, most of the beads will be deep within the rigid |
475 |
+ |
body and will not contribute to the hydrodynamic tensor. In the {\it |
476 |
+ |
rough shell} approach, beads which have all of their lattice neighbors |
477 |
+ |
inside the structure are considered interior beads, and are removed |
478 |
+ |
from the calculation. After following this procedure, only those |
479 |
+ |
beads in direct contact with the van der Waals surface of the rigid |
480 |
+ |
body are retained. For reasonably large molecular structures, this |
481 |
+ |
truncation can still produce bead assemblies with thousands of |
482 |
+ |
members. |
483 |
+ |
|
484 |
+ |
If all of the {\it atoms} comprising the rigid substructure are |
485 |
+ |
spherical and non-overlapping, the tensor in |
486 |
+ |
Eq.~(\ref{introEquation:RPTensorNonOverlapped}) may be used directly |
487 |
+ |
using the atoms themselves as the hydrodynamic beads. This is a |
488 |
+ |
variant of the {\it bead model} approach of Carrasco and Garc\'{i}a de |
489 |
+ |
la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$ |
490 |
+ |
matrix can be quite small, and the calculation of the hydrodynamic |
491 |
+ |
tensor is straightforward. |
492 |
+ |
|
493 |
+ |
In general, the inversion of the ${\bf B}$ matrix is the most |
494 |
+ |
computationally demanding task. This inversion is done only once for |
495 |
+ |
each type of rigid structure. We have used straightforward |
496 |
+ |
LU-decomposition to solve the linear system and to obtain the elements |
497 |
+ |
of ${\bf C}$. Once ${\bf C}$ has been obtained, the location of the |
498 |
+ |
center of resistance ($R$) is found and the resistance tensor at this |
499 |
+ |
point is calculated. The $3 \times 1$ vector giving the location of |
500 |
+ |
the rigid body's center of resistance and the $6 \times 6$ resistance |
501 |
+ |
tensor are then stored for use in the Langevin dynamics calculation. |
502 |
+ |
These quantities depend on solvent viscosity and temperature and must |
503 |
+ |
be recomputed if different simulation conditions are required. |
504 |
+ |
|
505 |
|
\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
506 |
|
|
507 |
|
Consider the Langevin equations of motion in generalized coordinates |
508 |
|
\begin{equation} |
509 |
< |
\mathbf{M} \dot{\mathbf{V}}(t) = \mathbf{F}_{s}(t) + |
510 |
< |
\mathbf{F}_{f}(t) + \mathbf{F}_{r}(t) |
509 |
> |
{\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) + |
510 |
> |
{\bf F}_{f}(t) + {\bf F}_{r}(t) |
511 |
|
\label{LDGeneralizedForm} |
512 |
|
\end{equation} |
513 |
< |
where $\mathbf{M}$ is a $6 \times 6$ diagonal mass matrix (which |
513 |
> |
where ${\bf M}$ is a $6 \times 6$ diagonal mass matrix (which |
514 |
|
includes the mass of the rigid body as well as the moments of inertia |
515 |
< |
in the body-fixed frame) and $\mathbf{V}$ is a generalized velocity, |
516 |
< |
$\mathbf{V} = |
517 |
< |
\left\{\mathbf{v},\mathbf{\omega}\right\}$. The right side of |
515 |
> |
in the body-fixed frame) and ${\bf V}$ is a generalized velocity, |
516 |
> |
${\bf V} = |
517 |
> |
\left\{{\bf v},{\bf \omega}\right\}$. The right side of |
518 |
|
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a |
519 |
< |
system force $\mathbf{F}_{s}$, a frictional or dissipative force |
520 |
< |
$\mathbf{F}_{f}$ and stochastic force $\mathbf{F}_{r}$. While the |
521 |
< |
evolution of the system in Newtonian mechanics is typically done in |
522 |
< |
the lab-fixed frame, it is convenient to handle the dynamics of rigid |
523 |
< |
bodies in the body-fixed frame. Thus the friction and random forces |
524 |
< |
are calculated in body-fixed frame and may be converted back to |
525 |
< |
lab-fixed frame using the rigid body's rotation matrix ($Q$): |
519 |
> |
system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf |
520 |
> |
F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution |
521 |
> |
of the system in Newtonian mechanics is typically done in the lab |
522 |
> |
frame, it is convenient to handle the dynamics of rigid bodies in |
523 |
> |
body-fixed frames. Thus the friction and random forces on each |
524 |
> |
substructure are calculated in a body-fixed frame and may converted |
525 |
> |
back to the lab frame using that substructure's rotation matrix (${\bf |
526 |
> |
Q}$): |
527 |
|
\begin{equation} |
528 |
< |
\mathbf{F}_{f,r} = |
528 |
> |
{\bf F}_{f,r} = |
529 |
|
\left( \begin{array}{c} |
530 |
< |
\mathbf{f}_{f,r} \\ |
531 |
< |
\mathbf{\tau}_{f,r} |
530 |
> |
{\bf f}_{f,r} \\ |
531 |
> |
{\bf \tau}_{f,r} |
532 |
|
\end{array} \right) |
533 |
|
= |
534 |
|
\left( \begin{array}{c} |
535 |
< |
Q^{T} \mathbf{f}^{b}_{f,r} \\ |
536 |
< |
Q^{T} \mathbf{\tau}^{b}_{f,r} |
535 |
> |
{\bf Q}^{T} {\bf f}^{~b}_{f,r} \\ |
536 |
> |
{\bf Q}^{T} {\bf \tau}^{~b}_{f,r} |
537 |
|
\end{array} \right) |
538 |
|
\end{equation} |
539 |
< |
The body-fixed friction force, $\mathbf{F}_{f}^b$, is proportional to |
540 |
< |
the velocity at the center of resistance $\mathbf{v}_{R}^b$ (in the |
541 |
< |
body-fixed frame) and the angular velocity $\mathbf{\omega}$ |
539 |
> |
The body-fixed friction force, ${\bf F}_{f}^{~b}$, is proportional to |
540 |
> |
the (body-fixed) velocity at the center of resistance |
541 |
> |
${\bf v}_{R}^{~b}$ and the angular velocity ${\bf \omega}$ |
542 |
|
\begin{equation} |
543 |
< |
\mathbf{F}_{f}^b (t) = \left( \begin{array}{l} |
544 |
< |
\mathbf{f}_{f}^b (t) \\ |
545 |
< |
\mathbf{\tau}_{f}^b (t) \\ |
543 |
> |
{\bf F}_{f}^{~b}(t) = \left( \begin{array}{l} |
544 |
> |
{\bf f}_{f}^{~b}(t) \\ |
545 |
> |
{\bf \tau}_{f}^{~b}(t) \\ |
546 |
|
\end{array} \right) = - \left( \begin{array}{*{20}c} |
547 |
< |
\Xi_{R,t} & \Xi_{R,c}^T \\ |
548 |
< |
\Xi_{R,c} & \Xi_{R,r} \\ |
547 |
> |
\Xi_{R}^{tt} & \Xi_{R}^{rt} \\ |
548 |
> |
\Xi_{R}^{tr} & \Xi_{R}^{rr} \\ |
549 |
|
\end{array} \right)\left( \begin{array}{l} |
550 |
< |
\mathbf{v}_{R}^b (t) \\ |
551 |
< |
\mathbf{\omega} (t) \\ |
550 |
> |
{\bf v}_{R}^{~b}(t) \\ |
551 |
> |
{\bf \omega}(t) \\ |
552 |
|
\end{array} \right), |
553 |
|
\end{equation} |
554 |
< |
while the random force, $\mathbf{F}_{r}$, is a Gaussian stochastic |
555 |
< |
variable with zero mean and variance |
554 |
> |
while the random force, ${\bf F}_{r}$, is a Gaussian stochastic |
555 |
> |
variable with zero mean and variance, |
556 |
|
\begin{equation} |
557 |
< |
\left\langle {\mathbf{F}_{r}(t) (\mathbf{F}_{r}(t'))^T } \right\rangle = |
558 |
< |
\left\langle {\mathbf{F}_{r}^b (t) (\mathbf{F}_{r}^b (t'))^T } \right\rangle = |
559 |
< |
2 k_B T \Xi_R \delta(t - t'). \label{randomForce} |
557 |
> |
\left\langle {{\bf F}_{r}(t) ({\bf F}_{r}(t'))^T } \right\rangle = |
558 |
> |
\left\langle {{\bf F}_{r}^{~b} (t) ({\bf F}_{r}^{~b} (t'))^T } \right\rangle = |
559 |
> |
2 k_B T \Xi_R \delta(t - t'). \label{eq:randomForce} |
560 |
|
\end{equation} |
561 |
|
$\Xi_R$ is the $6\times6$ resistance tensor at the center of |
562 |
< |
resistance. Once this tensor is known for a given rigid body, |
563 |
< |
obtaining a stochastic vector that has the properties in |
564 |
< |
Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a |
565 |
< |
one-time Cholesky decomposition to obtain the square root matrix of |
566 |
< |
the resistance tensor $\Xi_R = \mathbf{S} \mathbf{S}^{T}$, where |
567 |
< |
$\mathbf{S}$ is a lower triangular matrix.\cite{SchlickBook} A vector |
568 |
< |
with the statistics required for the random force can then be obtained |
569 |
< |
by multiplying $\mathbf{S}$ onto a 6-vector $Z$ which has elements |
570 |
< |
chosen from a Gaussian distribution, such that: |
562 |
> |
resistance. Once this tensor is known for a given rigid body (as |
563 |
> |
described in the previous section) obtaining a stochastic vector that |
564 |
> |
has the properties in Eq. (\ref{eq:randomForce}) can be done |
565 |
> |
efficiently by carrying out a one-time Cholesky decomposition to |
566 |
> |
obtain the square root matrix of the resistance tensor, |
567 |
> |
\begin{equation} |
568 |
> |
\Xi_R = {\bf S} {\bf S}^{T}, |
569 |
> |
\label{eq:Cholesky} |
570 |
> |
\end{equation} |
571 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
572 |
> |
vector with the statistics required for the random force can then be |
573 |
> |
obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which |
574 |
> |
has elements chosen from a Gaussian distribution, such that: |
575 |
|
\begin{equation} |
576 |
< |
\langle Z_i \rangle = 0, \hspace{1in} \langle Z_i \cdot Z_j \rangle = \frac{2 k_B |
577 |
< |
T}{\delta t} \delta_{ij}. |
576 |
> |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
577 |
> |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
578 |
|
\end{equation} |
579 |
< |
The random force, $F_{r}^{b} = \mathbf{S} Z$, can be shown to have the |
580 |
< |
correct ohmic |
579 |
> |
where $\delta t$ is the timestep in use during the simulation. The |
580 |
> |
random force, ${\bf F}_{r}^{~b} = {\bf S} {\bf Z}$, can be shown to have the |
581 |
> |
correct properties required by Eq. (\ref{eq:randomForce}). |
582 |
|
|
583 |
< |
|
584 |
< |
Each |
530 |
< |
time a random force vector is needed, a gaussian random vector is |
531 |
< |
generated and then the square root matrix is multiplied onto this |
532 |
< |
vector. |
533 |
< |
|
534 |
< |
The equation of motion for $\mathbf{v}$ can be written as |
583 |
> |
The equation of motion for the translational velocity of the center of |
584 |
> |
mass (${\bf v}$) can be written as |
585 |
|
\begin{equation} |
586 |
< |
m \dot{\mathbf{v}} (t) = \mathbf{f}_{s}^l (t) + \mathbf{f}_{f}^l (t) + |
587 |
< |
\mathbf{f}_{r}^l (t) |
586 |
> |
m \dot{{\bf v}} (t) = {\bf f}_{s}(t) + {\bf f}_{f}(t) + |
587 |
> |
{\bf f}_{r}(t) |
588 |
|
\end{equation} |
589 |
< |
Since the frictional force is applied at the center of resistance |
590 |
< |
which generally does not coincide with the center of mass, an extra |
591 |
< |
torque is exerted at the center of mass. Thus, the net body-fixed |
592 |
< |
frictional torque at the center of mass, $\tau_{f}^b (t)$, is |
593 |
< |
given by |
589 |
> |
Since the frictional and random forces are applied at the center of |
590 |
> |
resistance, which generally does not coincide with the center of mass, |
591 |
> |
extra torques are exerted at the center of mass. Thus, the net |
592 |
> |
body-fixed torque at the center of mass, $\tau^{~b}(t)$, |
593 |
> |
is given by |
594 |
|
\begin{equation} |
595 |
< |
\tau_{f}^b \leftarrow \tau_{f}^b + \mathbf{r}_{MR} \times \mathbf{f}_{f}^b |
595 |
> |
\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) |
596 |
|
\end{equation} |
597 |
< |
where $r_{MR}$ is the vector from the center of mass to the center |
598 |
< |
of the resistance. Instead of integrating the angular velocity in |
599 |
< |
lab-fixed frame, we consider the equation of angular momentum in |
600 |
< |
body-fixed frame |
597 |
> |
where ${\bf r}_{MR}$ is the vector from the center of mass to the center of |
598 |
> |
resistance. Instead of integrating the angular velocity in lab-fixed |
599 |
> |
frame, we consider the equation of motion for the angular momentum |
600 |
> |
(${\bf j}$) in the body-fixed frame |
601 |
|
\begin{equation} |
602 |
< |
\dot j(t) = \tau_{s}^b (t) + \tau_{f}^b (t) + \tau_{r}^b(t) |
602 |
> |
\frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t) |
603 |
|
\end{equation} |
604 |
< |
Embedding the friction terms into force and torque, one can integrate |
605 |
< |
the Langevin equations of motion for rigid body of arbitrary shape in |
606 |
< |
a velocity-Verlet style 2-part algorithm, where $h= \delta t$: |
604 |
> |
Embedding the friction and random forces into the the total force and |
605 |
> |
torque, one can integrate the Langevin equations of motion for a rigid |
606 |
> |
body of arbitrary shape in a velocity-Verlet style 2-part algorithm, |
607 |
> |
where $h = \delta t$: |
608 |
|
|
609 |
< |
{\tt moveA:} |
609 |
> |
{\tt move A:} |
610 |
|
\begin{align*} |
611 |
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
612 |
|
+ \frac{h}{2} \left( {\bf f}(t) / m \right), \\ |
615 |
|
+ h {\bf v}\left(t + h / 2 \right), \\ |
616 |
|
% |
617 |
|
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
618 |
< |
+ \frac{h}{2} {\bf \tau}^b(t), \\ |
618 |
> |
+ \frac{h}{2} {\bf \tau}^{~b}(t), \\ |
619 |
|
% |
620 |
< |
\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
620 |
> |
{\bf Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
621 |
|
(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right). |
622 |
|
\end{align*} |
623 |
< |
In this context, the $\mathrm{rotate}$ function is the reversible |
624 |
< |
product of the three body-fixed rotations, |
623 |
> |
In this context, $\overleftrightarrow{\mathsf{I}}$ is the diagonal |
624 |
> |
moment of inertia tensor, and the $\mathrm{rotate}$ function is the |
625 |
> |
reversible product of the three body-fixed rotations, |
626 |
|
\begin{equation} |
627 |
|
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
628 |
|
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y |
629 |
|
/ 2) \cdot \mathsf{G}_x(a_x /2), |
630 |
|
\end{equation} |
631 |
|
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, |
632 |
< |
rotates both the rotation matrix ($\mathsf{Q}$) and the body-fixed |
632 |
> |
rotates both the rotation matrix ($\mathbf{Q}$) and the body-fixed |
633 |
|
angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed |
634 |
|
axis $\alpha$, |
635 |
|
\begin{equation} |
636 |
|
\mathsf{G}_\alpha( \theta ) = \left\{ |
637 |
|
\begin{array}{lcl} |
638 |
< |
\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
638 |
> |
\mathbf{Q}(t) & \leftarrow & \mathbf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
639 |
|
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf |
640 |
|
j}(0). |
641 |
|
\end{array} |
657 |
|
\end{equation} |
658 |
|
All other rotations follow in a straightforward manner. After the |
659 |
|
first part of the propagation, the forces and body-fixed torques are |
660 |
< |
calculated at the new positions and orientations |
660 |
> |
calculated at the new positions and orientations. The system forces |
661 |
> |
and torques are derivatives of the total potential energy function |
662 |
> |
($U$) with respect to the rigid body positions (${\bf r}$) and the |
663 |
> |
columns of the transposed rotation matrix ${\bf Q}^T = \left({\bf |
664 |
> |
u}_x, {\bf u}_y, {\bf u}_z \right)$: |
665 |
|
|
666 |
< |
{\tt doForces:} |
666 |
> |
{\tt Forces:} |
667 |
|
\begin{align*} |
668 |
< |
{\bf f}(t + h) &\leftarrow |
669 |
< |
- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\ |
668 |
> |
{\bf f}_{s}(t + h) & \leftarrow |
669 |
> |
- \left(\frac{\partial U}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\ |
670 |
|
% |
671 |
< |
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) |
672 |
< |
\times \frac{\partial V}{\partial {\bf u}}, \\ |
671 |
> |
{\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h) |
672 |
> |
\times \frac{\partial U}{\partial {\bf u}} \\ |
673 |
|
% |
674 |
< |
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h) |
675 |
< |
\cdot {\bf \tau}^s(t + h). |
674 |
> |
{\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) \\ |
675 |
> |
% |
676 |
> |
{\bf f}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tt} \cdot |
677 |
> |
{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rt} \cdot {\bf \omega}(t+h) \\ |
678 |
> |
% |
679 |
> |
{\bf \tau}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tr} \cdot |
680 |
> |
{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rr} \cdot {\bf \omega}(t+h) \\ |
681 |
> |
% |
682 |
> |
Z & \leftarrow {\tt GaussianNormal}(2 k_B T / h, 6) \\ |
683 |
> |
{\bf F}_{R,r}^{b}(t+h) & \leftarrow {\bf S} \cdot Z \\ |
684 |
> |
% |
685 |
> |
{\bf f}(t+h) & \leftarrow {\bf f}_{s}(t+h) + \mathbf{Q}^{T}(t+h) |
686 |
> |
\cdot \left({\bf f}_{R,f}^{~b} + {\bf f}_{R,r}^{~b} \right) \\ |
687 |
> |
% |
688 |
> |
\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) \\ |
689 |
> |
\tau^{~b}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \tau(t+h) \\ |
690 |
|
\end{align*} |
691 |
+ |
Frictional (and random) forces and torques must be computed at the |
692 |
+ |
center of resistance, so there are additional steps required to find |
693 |
+ |
the body-fixed velocity (${\bf v}_{R}^{~b}$) at this location. Mapping |
694 |
+ |
the frictional and random forces at the center of resistance back to |
695 |
+ |
the center of mass also introduces an additional term in the torque |
696 |
+ |
one obtains at the center of mass. |
697 |
+ |
|
698 |
|
Once the forces and torques have been obtained at the new time step, |
699 |
|
the velocities can be advanced to the same time value. |
700 |
|
|
701 |
< |
{\tt moveB:} |
701 |
> |
{\tt move B:} |
702 |
|
\begin{align*} |
703 |
|
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2 |
704 |
|
\right) |
706 |
|
% |
707 |
|
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2 |
708 |
|
\right) |
709 |
< |
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
709 |
> |
+ \frac{h}{2} {\bf \tau}^{~b}(t + h) . |
710 |
|
\end{align*} |
711 |
|
|
712 |
|
\section{Validating the Method\label{sec:validating}} |
769 |
|
We performed reference microcanonical simulations with explicit |
770 |
|
solvents for each of the different model system. In each case there |
771 |
|
was one solute model and 1929 solvent molecules present in the |
772 |
< |
simulation box. All simulations were equilibrated using a |
772 |
> |
simulation box. All simulations were equilibrated for 5 ns using a |
773 |
|
constant-pressure and temperature integrator with target values of 300 |
774 |
|
K for the temperature and 1 atm for pressure. Following this stage, |
775 |
< |
further equilibration and sampling was done in a microcanonical |
776 |
< |
ensemble. Since the model bodies are typically quite massive, we were |
777 |
< |
able to use a time step of 25 fs. |
775 |
> |
further equilibration (5 ns) and sampling (10 ns) was done in a |
776 |
> |
microcanonical ensemble. Since the model bodies are typically quite |
777 |
> |
massive, we were able to use a time step of 25 fs. |
778 |
|
|
779 |
|
The model systems studied used both Lennard-Jones spheres as well as |
780 |
|
uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne |
781 |
|
potential was a single site model for the interactions of rigid |
782 |
< |
ellipsoidal molecules.\cite{Gay81} It can be thought of as a |
782 |
> |
ellipsoidal molecules.\cite{Gay1981} It can be thought of as a |
783 |
|
modification of the Gaussian overlap model originally described by |
784 |
|
Berne and Pechukas.\cite{Berne72} The potential is constructed in the |
785 |
|
familiar form of the Lennard-Jones function using |
786 |
|
orientation-dependent $\sigma$ and $\epsilon$ parameters, |
787 |
|
\begin{equation*} |
788 |
< |
V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat |
789 |
< |
r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, |
790 |
< |
{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u |
788 |
> |
V_{ij}({{\bf \hat u}_i}, {{\bf \hat u}_j}, {{\bf \hat |
789 |
> |
r}_{ij}}) = 4\epsilon ({{\bf \hat u}_i}, {{\bf \hat u}_j}, |
790 |
> |
{{\bf \hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u |
791 |
|
}_i}, |
792 |
< |
{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12} |
793 |
< |
-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, |
794 |
< |
{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right] |
792 |
> |
{{\bf \hat u}_j}, {{\bf \hat r}_{ij}})+\sigma_0}\right)^{12} |
793 |
> |
-\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u}_i}, {{\bf \hat u}_j}, |
794 |
> |
{{\bf \hat r}_{ij}})+\sigma_0}\right)^6\right] |
795 |
|
\label{eq:gb} |
796 |
|
\end{equation*} |
797 |
|
|
808 |
|
Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e / |
809 |
|
\epsilon^s$, describes the ratio between the well depths in the {\it |
810 |
|
end-to-end} and side-by-side configurations. Details of the potential |
811 |
< |
are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an |
811 |
> |
are given elsewhere,\cite{Luckhurst90,Golubkov06,SunX._jp0762020} and an |
812 |
|
excellent overview of the computational methods that can be used to |
813 |
|
efficiently compute forces and torques for this potential can be found |
814 |
|
in Ref. \citen{Golubkov06} |
843 |
|
\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}. |
844 |
|
\label{eq:shear} |
845 |
|
\end{equation} |
846 |
< |
A similar form exists for the bulk viscosity |
847 |
< |
\begin{equation} |
771 |
< |
\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left( |
772 |
< |
\int_{t_0}^{t_0 + t} |
773 |
< |
\left(P\left(t'\right)-\left\langle P \right\rangle \right)dt' |
774 |
< |
\right)^2 \right\rangle_{t_0}. |
775 |
< |
\end{equation} |
776 |
< |
Alternatively, the shear viscosity can also be calculated using a |
777 |
< |
Green-Kubo formula with the off-diagonal pressure tensor correlation function, |
778 |
< |
\begin{equation} |
779 |
< |
\eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0 |
780 |
< |
+ t) \right\rangle_{t_0} dt, |
781 |
< |
\end{equation} |
782 |
< |
although this method converges extremely slowly and is not practical |
783 |
< |
for obtaining viscosities from molecular dynamics simulations. |
846 |
> |
which converges much more rapidly in molecular dynamics simulations |
847 |
> |
than the traditional Green-Kubo formula. |
848 |
|
|
849 |
|
The Langevin dynamics for the different model systems were performed |
850 |
|
at the same temperature as the average temperature of the |
870 |
|
compute the diffusive behavior for motion parallel to each body-fixed |
871 |
|
axis by projecting the displacement of the particle onto the |
872 |
|
body-fixed reference frame at $t=0$. With an isotropic solvent, as we |
873 |
< |
have used in this study, there are differences between the three |
874 |
< |
diffusion constants, but these must converge to the same value at |
875 |
< |
longer times. Translational diffusion constants for the different |
876 |
< |
shaped models are shown in table \ref{tab:translation}. |
873 |
> |
have used in this study, there may be differences between the three |
874 |
> |
diffusion constants at short times, but these must converge to the |
875 |
> |
same value at longer times. Translational diffusion constants for the |
876 |
> |
different shaped models are shown in table \ref{tab:translation}. |
877 |
|
|
878 |
|
In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational |
879 |
|
diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object |
953 |
|
an arbitrary value of 0.8 kcal/mol. |
954 |
|
|
955 |
|
The Stokes-Einstein behavior of large spherical particles in |
956 |
< |
hydrodynamic flows is well known, giving translational friction |
957 |
< |
coefficients of $6 \pi \eta R$ (stick boundary conditions) and |
958 |
< |
rotational friction coefficients of $8 \pi \eta R^3$. Recently, |
959 |
< |
Schmidt and Skinner have computed the behavior of spherical tag |
960 |
< |
particles in molecular dynamics simulations, and have shown that {\it |
961 |
< |
slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more |
962 |
< |
appropriate for molecule-sized spheres embedded in a sea of spherical |
963 |
< |
solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
956 |
> |
hydrodynamic flows with ``stick'' boundary conditions is well known, |
957 |
> |
and is given in Eqs. (\ref{eq:StokesTranslation}) and |
958 |
> |
(\ref{eq:StokesRotation}). Recently, Schmidt and Skinner have |
959 |
> |
computed the behavior of spherical tag particles in molecular dynamics |
960 |
> |
simulations, and have shown that {\it slip} boundary conditions |
961 |
> |
($\Xi_{tt} = 4 \pi \eta \rho$) may be more appropriate for |
962 |
> |
molecule-sized spheres embedded in a sea of spherical solvent |
963 |
> |
particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
964 |
|
|
965 |
|
Our simulation results show similar behavior to the behavior observed |
966 |
|
by Schmidt and Skinner. The diffusion constant obtained from our |
985 |
|
can be combined to give a single translational diffusion |
986 |
|
constant,\cite{Berne90} |
987 |
|
\begin{equation} |
988 |
< |
D = \frac{k_B T}{6 \pi \eta a} G(\rho), |
988 |
> |
D = \frac{k_B T}{6 \pi \eta a} G(s), |
989 |
|
\label{Dperrin} |
990 |
|
\end{equation} |
991 |
|
as well as a single rotational diffusion coefficient, |
992 |
|
\begin{equation} |
993 |
< |
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2) |
994 |
< |
G(\rho) - 1}{1 - \rho^4} \right\}. |
993 |
> |
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - s^2) |
994 |
> |
G(s) - 1}{1 - s^4} \right\}. |
995 |
|
\label{ThetaPerrin} |
996 |
|
\end{equation} |
997 |
< |
In these expressions, $G(\rho)$ is a function of the axial ratio |
998 |
< |
($\rho = b / a$), which for prolate ellipsoids, is |
997 |
> |
In these expressions, $G(s)$ is a function of the axial ratio |
998 |
> |
($s = b / a$), which for prolate ellipsoids, is |
999 |
|
\begin{equation} |
1000 |
< |
G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 - |
937 |
< |
\rho^2)^{1/2}}{\rho} \right\} |
1000 |
> |
G(s) = (1- s^2)^{-1/2} \ln \left\{ \frac{1 + (1 - s^2)^{1/2}}{s} \right\} |
1001 |
|
\label{GPerrin} |
1002 |
|
\end{equation} |
1003 |
|
Again, there is some uncertainty about the correct boundary conditions |
1022 |
|
exact treatment of the diffusion tensor as well as the rough-shell |
1023 |
|
model for the ellipsoid. |
1024 |
|
|
1025 |
< |
The translational diffusion constants from the microcanonical simulations |
1026 |
< |
agree well with the predictions of the Perrin model, although the rotational |
1027 |
< |
correlation times are a factor of 2 shorter than expected from hydrodynamic |
1028 |
< |
theory. One explanation for the slower rotation |
1029 |
< |
of explicitly-solvated ellipsoids is the possibility that solute-solvent |
1030 |
< |
collisions happen at both ends of the solute whenever the principal |
1031 |
< |
axis of the ellipsoid is turning. In the upper portion of figure |
1032 |
< |
\ref{fig:explanation} we sketch a physical picture of this explanation. |
1033 |
< |
Since our Langevin integrator is providing nearly quantitative agreement with |
1034 |
< |
the Perrin model, it also predicts orientational diffusion for ellipsoids that |
1035 |
< |
exceed explicitly solvated correlation times by a factor of two. |
1025 |
> |
The translational diffusion constants from the microcanonical |
1026 |
> |
simulations agree well with the predictions of the Perrin model, |
1027 |
> |
although the {\it rotational} correlation times are a factor of 2 |
1028 |
> |
shorter than expected from hydrodynamic theory. One explanation for |
1029 |
> |
the slower rotation of explicitly-solvated ellipsoids is the |
1030 |
> |
possibility that solute-solvent collisions happen at both ends of the |
1031 |
> |
solute whenever the principal axis of the ellipsoid is turning. In the |
1032 |
> |
upper portion of figure \ref{fig:explanation} we sketch a physical |
1033 |
> |
picture of this explanation. Since our Langevin integrator is |
1034 |
> |
providing nearly quantitative agreement with the Perrin model, it also |
1035 |
> |
predicts orientational diffusion for ellipsoids that exceed explicitly |
1036 |
> |
solvated correlation times by a factor of two. |
1037 |
|
|
1038 |
|
\subsection{Rigid dumbbells} |
1039 |
|
Perhaps the only {\it composite} rigid body for which analytic |
1041 |
|
two-sphere dumbbell model. This model consists of two non-overlapping |
1042 |
|
spheres held by a rigid bond connecting their centers. There are |
1043 |
|
competing expressions for the 6x6 resistance tensor for this |
1044 |
< |
model. Equation (\ref{introEquation:oseenTensor}) above gives the |
1045 |
< |
original Oseen tensor, while the second order expression introduced by |
1046 |
< |
Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la |
983 |
< |
Torre and Bloomfield,\cite{Torre1977} is given above as |
1044 |
> |
model. The second order expression introduced by Rotne and |
1045 |
> |
Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la Torre and |
1046 |
> |
Bloomfield,\cite{Torre1977} is given above as |
1047 |
|
Eq. (\ref{introEquation:RPTensorNonOverlapped}). In our case, we use |
1048 |
|
a model dumbbell in which the two spheres are identical Lennard-Jones |
1049 |
|
particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at |
1056 |
|
motion in a flow {\it perpendicular} to the inter-sphere |
1057 |
|
axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it |
1058 |
|
orientational} correlation times for this model system (other than |
1059 |
< |
those derived from the 6 x 6 tensors mentioned above). |
1059 |
> |
those derived from the 6 x 6 tensor mentioned above). |
1060 |
|
|
1061 |
|
The bead model for this model system comprises the two large spheres |
1062 |
|
by themselves, while the rough shell approximation used 3368 separate |
1124 |
|
} \label{fig:explanation} |
1125 |
|
\end{figure} |
1126 |
|
|
1064 |
– |
|
1065 |
– |
|
1127 |
|
\subsection{Composite banana-shaped molecules} |
1128 |
|
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids have |
1129 |
|
been used by Orlandi {\it et al.} to observe mesophases in |
1135 |
|
behavior of this model, we have left out the dipolar interactions of |
1136 |
|
the original Orlandi model. |
1137 |
|
|
1138 |
< |
A reference system composed of a single banana rigid body embedded in a |
1139 |
< |
sea of 1929 solvent particles was created and run under standard |
1140 |
< |
(microcanonical) molecular dynamics. The resulting viscosity of this |
1141 |
< |
mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})). |
1142 |
< |
To calculate the hydrodynamic properties of the banana rigid body model, |
1143 |
< |
we created a rough shell (see Fig.~\ref{fig:roughShell}), in which |
1144 |
< |
the banana is represented as a ``shell'' made of 3321 identical beads |
1145 |
< |
(0.25 \AA\ in diameter) distributed on the surface. Applying the |
1146 |
< |
procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
1147 |
< |
identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 \AA), as |
1148 |
< |
well as the resistance tensor, |
1149 |
< |
\begin{equation*} |
1089 |
< |
\Xi = |
1090 |
< |
\left( {\begin{array}{*{20}c} |
1091 |
< |
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
1092 |
< |
0& 0.9270&-0.007063& 0.08585&0&0\\ |
1093 |
< |
0&-0.007063&0.7494&0.2057&0&0\\ |
1094 |
< |
0&0.0858&0.2057& 58.64& 0&0\\0.08585&0&0&0&48.30&3.219&\\0.2057&0&0&0&3.219&10.7373\\\end{array}} \right), |
1095 |
< |
\end{equation*} |
1096 |
< |
where the units for translational, translation-rotation coupling and |
1097 |
< |
rotational tensors are (kcal fs / mol \AA$^2$), (kcal fs / mol \AA\ rad), |
1098 |
< |
and (kcal fs / mol rad$^2$), respectively. |
1138 |
> |
A reference system composed of a single banana rigid body embedded in |
1139 |
> |
a sea of 1929 solvent particles was created and run under standard |
1140 |
> |
(microcanonical) molecular dynamics. The resulting viscosity of this |
1141 |
> |
mixture was 0.298 centipoise (as estimated using |
1142 |
> |
Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties of |
1143 |
> |
the banana rigid body model, we created a rough shell (see |
1144 |
> |
Fig.~\ref{fig:roughShell}), in which the banana is represented as a |
1145 |
> |
``shell'' made of 3321 identical beads (0.25 \AA\ in diameter) |
1146 |
> |
distributed on the surface. Applying the procedure described in |
1147 |
> |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
1148 |
> |
identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 |
1149 |
> |
\AA). |
1150 |
|
|
1151 |
< |
The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor) |
1152 |
< |
are essentially quantitative for translational diffusion of this model. |
1153 |
< |
Orientational correlation times under the Langevin rigid-body integrator |
1154 |
< |
are within 11\% of the values obtained from explicit solvent, but these |
1155 |
< |
models also exhibit some solvent inaccessible surface area in the |
1156 |
< |
explicitly-solvated case. |
1151 |
> |
The Langevin rigid-body integrator (and the hydrodynamic diffusion |
1152 |
> |
tensor) are essentially quantitative for translational diffusion of |
1153 |
> |
this model. Orientational correlation times under the Langevin |
1154 |
> |
rigid-body integrator are within 11\% of the values obtained from |
1155 |
> |
explicit solvent, but these models also exhibit some solvent |
1156 |
> |
inaccessible surface area in the explicitly-solvated case. |
1157 |
|
|
1158 |
|
\subsection{Composite sphero-ellipsoids} |
1159 |
+ |
|
1160 |
|
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
1161 |
|
used recently as models for lipid |
1162 |
< |
molecules.\cite{SunGezelter08,Ayton01} |
1163 |
< |
MORE DETAILS |
1164 |
< |
|
1165 |
< |
A reference system composed of a single lipid rigid body embedded in a |
1114 |
< |
sea of 1929 solvent particles was created and run under standard |
1115 |
< |
(microcanonical) molecular dynamics. The resulting viscosity of this |
1116 |
< |
mixture was 0.349 centipoise (as estimated using |
1162 |
> |
molecules.\cite{SunX._jp0762020,Ayton01} A reference system composed of |
1163 |
> |
a single lipid rigid body embedded in a sea of 1929 solvent particles |
1164 |
> |
was created and run under a microcanonical ensemble. The resulting |
1165 |
> |
viscosity of this mixture was 0.349 centipoise (as estimated using |
1166 |
|
Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties of |
1167 |
|
the lipid rigid body model, we created a rough shell (see |
1168 |
|
Fig.~\ref{fig:roughShell}), in which the lipid is represented as a |
1172 |
|
identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46 |
1173 |
|
\AA). |
1174 |
|
|
1175 |
+ |
The translational diffusion constants and rotational correlation times |
1176 |
+ |
obtained using the Langevin rigid-body integrator (and the |
1177 |
+ |
hydrodynamic tensor) are essentially quantitative when compared with |
1178 |
+ |
the explicit solvent simulations for this model system. |
1179 |
|
|
1180 |
< |
\subsection{Summary} |
1181 |
< |
According to our simulations, the langevin dynamics is a reliable |
1182 |
< |
theory to apply to replace the explicit solvents, especially for the |
1183 |
< |
translation properties. For large molecules, the rotation properties |
1184 |
< |
are also mimiced reasonablly well. |
1180 |
> |
\subsection{Summary of comparisons with explicit solvent simulations} |
1181 |
> |
The Langevin rigid-body integrator we have developed is a reliable way |
1182 |
> |
to replace explicit solvent simulations in cases where the detailed |
1183 |
> |
solute-solvent interactions do not greatly impact the behavior of the |
1184 |
> |
solute. As such, it has the potential to greatly increase the length |
1185 |
> |
and time scales of coarse grained simulations of large solvated |
1186 |
> |
molecules. In cases where the dielectric screening of the solvent, or |
1187 |
> |
specific solute-solvent interactions become important for structural |
1188 |
> |
or dynamic features of the solute molecule, this integrator may be |
1189 |
> |
less useful. However, for the kinds of coarse-grained modeling that |
1190 |
> |
have become popular in recent years (ellipsoidal side chains, rigid |
1191 |
> |
bodies, and molecular-scale models), this integrator may prove itself |
1192 |
> |
to be quite valuable. |
1193 |
|
|
1194 |
|
\begin{figure} |
1195 |
|
\centering |
1215 |
|
\caption{Translational diffusion constants (D) for the model systems |
1216 |
|
calculated using microcanonical simulations (with explicit solvent), |
1217 |
|
theoretical predictions, and Langevin simulations (with implicit solvent). |
1218 |
< |
Analytical solutions for the exactly-solved hydrodynamics models are |
1219 |
< |
from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936} |
1218 |
> |
Analytical solutions for the exactly-solved hydrodynamics models are obtained |
1219 |
> |
from: Stokes' law (sphere), and Refs. \citen{Perrin1934} and \citen{Perrin1936} |
1220 |
|
(ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq} |
1221 |
|
(dumbbell). The other model systems have no known analytic solution. |
1222 |
< |
All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= |
1222 |
> |
All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= |
1223 |
|
$10^{-4}$ \AA$^2$ / fs). } |
1224 |
|
\begin{tabular}{lccccccc} |
1225 |
|
\hline |
1233 |
|
dumbbell & 0.308 & 2.06 & & 1.64 & bead model & 1.65 & 1.62 \\ |
1234 |
|
& 0.308 & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\ |
1235 |
|
banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\ |
1236 |
< |
lipid & 0.349 & 0.96 & & & rough shell & 1.33 & 1.32 \\ |
1236 |
> |
lipid & 0.349 & 1.41 & & & rough shell & 1.33 & 1.32 \\ |
1237 |
|
\end{tabular} |
1238 |
|
\label{tab:translation} |
1239 |
|
\end{center} |
1272 |
|
|
1273 |
|
\section{Application: A rigid-body lipid bilayer} |
1274 |
|
|
1275 |
< |
The Langevin dynamics integrator was applied to study the formation of |
1276 |
< |
corrugated structures emerging from simulations of the coarse grained |
1277 |
< |
lipid molecular models presented above. The initial configuration is |
1278 |
< |
taken from our molecular dynamics studies on lipid bilayers with |
1279 |
< |
lennard-Jones sphere solvents. The solvent molecules were excluded |
1280 |
< |
from the system, and the experimental value for the viscosity of water |
1281 |
< |
at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects |
1282 |
< |
of the solvent. The absence of explicit solvent molecules and the |
1283 |
< |
stability of the integrator allowed us to take timesteps of 50 fs. A |
1284 |
< |
total simulation run time of 100 ns was sampled. |
1285 |
< |
Fig. \ref{fig:bilayer} shows the configuration of the system after 100 |
1286 |
< |
ns, and the ripple structure remains stable during the entire |
1287 |
< |
trajectory. Compared with using explicit bead-model solvent |
1288 |
< |
molecules, the efficiency of the simulation has increased by an order |
1289 |
< |
of magnitude. |
1275 |
> |
To test the accuracy and efficiency of the new integrator, we applied |
1276 |
> |
it to study the formation of corrugated structures emerging from |
1277 |
> |
simulations of the coarse grained lipid molecular models presented |
1278 |
> |
above. The initial configuration is taken from earlier molecular |
1279 |
> |
dynamics studies on lipid bilayers which had used spherical |
1280 |
> |
(Lennard-Jones) solvent particles and moderate (480 solvated lipid |
1281 |
> |
molecules) system sizes.\cite{SunX._jp0762020} the solvent molecules |
1282 |
> |
were excluded from the system and the box was replicated three times |
1283 |
> |
in the x- and y- axes (giving a single simulation cell containing 4320 |
1284 |
> |
lipids). The experimental value for the viscosity of water at 20C |
1285 |
> |
($\eta = 1.00$ cp) was used with the Langevin integrator to mimic the |
1286 |
> |
hydrodynamic effects of the solvent. The absence of explicit solvent |
1287 |
> |
molecules and the stability of the integrator allowed us to take |
1288 |
> |
timesteps of 50 fs. A simulation run time of 30 ns was sampled to |
1289 |
> |
calculate structural properties. Fig. \ref{fig:bilayer} shows the |
1290 |
> |
configuration of the system after 30 ns. Structural properties of the |
1291 |
> |
bilayer (e.g. the head and body $P_2$ order parameters) are nearly |
1292 |
> |
identical to those obtained via solvated molecular dynamics. The |
1293 |
> |
ripple structure remained stable during the entire trajectory. |
1294 |
> |
Compared with using explicit bead-model solvent molecules, the 30 ns |
1295 |
> |
trajectory for 4320 lipids with the Langevin integrator is now {\it |
1296 |
> |
faster} on the same hardware than the same length trajectory was for |
1297 |
> |
the 480-lipid system previously studied. |
1298 |
|
|
1299 |
|
\begin{figure} |
1300 |
|
\centering |
1301 |
|
\includegraphics[width=\linewidth]{bilayer} |
1302 |
|
\caption[Snapshot of a bilayer of rigid-body models for lipids]{A |
1303 |
< |
snapshot of a bilayer composed of rigid-body models for lipid |
1303 |
> |
snapshot of a bilayer composed of 4320 rigid-body models for lipid |
1304 |
|
molecules evolving using the Langevin integrator described in this |
1305 |
|
work.} \label{fig:bilayer} |
1306 |
|
\end{figure} |
1307 |
|
|
1308 |
|
\section{Conclusions} |
1309 |
|
|
1310 |
< |
We have presented a new Langevin algorithm by incorporating the |
1311 |
< |
hydrodynamics properties of arbitrary shaped molecules into an |
1312 |
< |
advanced symplectic integration scheme. Further studies in systems |
1313 |
< |
involving banana shaped molecules illustrated that the dynamic |
1314 |
< |
properties could be preserved by using this new algorithm as an |
1315 |
< |
implicit solvent model. |
1310 |
> |
We have presented a new algorithm for carrying out Langevin dynamics |
1311 |
> |
simulations on complex rigid bodies by incorporating the hydrodynamic |
1312 |
> |
resistance tensors for arbitrary shapes into an advanced symplectic |
1313 |
> |
integration scheme. The integrator gives quantitative agreement with |
1314 |
> |
both analytic and approximate hydrodynamic theories, and works |
1315 |
> |
reasonably well at reproducing the solute dynamical properties |
1316 |
> |
(diffusion constants, and orientational relaxation times) from |
1317 |
> |
explicitly-solvated simulations. For the cases where there are |
1318 |
> |
discrepancies between our Langevin integrator and the explicit solvent |
1319 |
> |
simulations, two features of molecular simulations help explain the |
1320 |
> |
differences. |
1321 |
|
|
1322 |
+ |
First, the use of ``stick'' boundary conditions for molecular-sized |
1323 |
+ |
solutes in a sea of similarly-sized solvent particles may be |
1324 |
+ |
problematic. We are certainly not the first group to notice this |
1325 |
+ |
difference between hydrodynamic theories and explicitly-solvated |
1326 |
+ |
molecular |
1327 |
+ |
simulations.\cite{Schmidt:2004fj,Schmidt:2003kx,Ravichandran:1999fk,TANG:1993lr} |
1328 |
+ |
The problem becomes particularly noticable in both the translational |
1329 |
+ |
diffusion of the spherical particles and the rotational diffusion of |
1330 |
+ |
the ellipsoids. In both of these cases it is clear that the |
1331 |
+ |
approximations that go into hydrodynamics are the source of the error, |
1332 |
+ |
and not the integrator itself. |
1333 |
|
|
1334 |
+ |
Second, in the case of structures which have substantial surface area |
1335 |
+ |
that is inaccessible to solvent particles, the hydrodynamic theories |
1336 |
+ |
(and the Langevin integrator) may overestimate the effects of solvent |
1337 |
+ |
friction because they overestimate the exposed surface area of the |
1338 |
+ |
rigid body. This is particularly noticable in the rotational |
1339 |
+ |
diffusion of the dumbbell model. We believe that given a solvent of |
1340 |
+ |
known radius, it may be possible to modify the rough shell approach to |
1341 |
+ |
place beads on solvent-accessible surface, instead of on the geometric |
1342 |
+ |
surface defined by the van der Waals radii of the components of the |
1343 |
+ |
rigid body. Further work to confirm the behavior of this new |
1344 |
+ |
approximation is ongoing. |
1345 |
+ |
|
1346 |
|
\section{Acknowledgments} |
1347 |
|
Support for this project was provided by the National Science |
1348 |
|
Foundation under grant CHE-0134881. T.L. also acknowledges the |
1349 |
< |
financial support from center of applied mathematics at University |
1350 |
< |
of Notre Dame. |
1349 |
> |
financial support from Center of Applied Mathematics at University of |
1350 |
> |
Notre Dame. |
1351 |
> |
|
1352 |
> |
\end{doublespace} |
1353 |
|
\newpage |
1354 |
|
|
1355 |
< |
\bibliographystyle{jcp} |
1355 |
> |
\bibliographystyle{jcp2} |
1356 |
|
\bibliography{langevin} |
1258 |
– |
|
1357 |
|
\end{document} |