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