| 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's theory, Klimov and Thirumalai identified free-energy |
| 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 |
| 85 |
|
finite structures.\cite{Berkov2005a}] |
| 86 |
|
|
| 87 |
|
In typical LD simulations, the friction and random forces on |
| 88 |
< |
individual atoms are taken from the Stokes-Einstein hydrodynamic |
| 89 |
< |
approximation, |
| 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 \\ |
| 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 |
< |
The use of rigid substructures,\cite{???} |
| 98 |
< |
coarse-graining,\cite{Ayton,Sun,Zannoni} and ellipsoidal |
| 99 |
< |
representations of protein side chains~\cite{Schulten} has made the |
| 100 |
< |
use of the Stokes-Einstein approximation problematic. A rigid |
| 101 |
< |
substructure moves as a single unit with orientational as well as |
| 102 |
< |
translational degrees of freedom. This requires a more general |
| 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{MarshallNewton,GarciaDeLaTorre} There a need to have a |
| 111 |
< |
more thorough treatment of hydrodynamics included in the library of |
| 112 |
< |
methods available for performing Langevin simulations. |
| 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 simulations.\cite{Sun2008} |
| 143 |
< |
Many of the water models in common use are also rigid-body |
| 144 |
< |
models,\cite{TIPs,SPC/E} although they are typically evolved using |
| 145 |
< |
constraints rather than rigid body equations of motion. |
| 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 |
| 150 |
< |
degrees of freedom. However, due to $1 \over \sin \theta$ |
| 151 |
< |
singularities, the numerical integration of corresponding equations of |
| 152 |
< |
these motion can become inaccurate (and inefficient). Although an |
| 153 |
< |
alternative integrator using multiple sets of Euler angles can |
| 154 |
< |
overcome this problem,\cite{Barojas1973} the computational penalty and |
| 155 |
< |
the loss of angular momentum conservation remain. A singularity-free |
| 156 |
< |
representation utilizing quaternions was developed by Evans in |
| 157 |
< |
1977.\cite{Evans1977} Unfortunately, this approach uses a nonseparable |
| 158 |
< |
Hamiltonian resulting from the quaternion representation, which |
| 159 |
< |
prevented symplectic algorithms from being utilized until very |
| 134 |
< |
recently.\cite{Miller2002} Another approach is the application of |
| 135 |
< |
holonomic constraints to the atoms belonging to the rigid body. Each |
| 136 |
< |
atom moves independently under the normal forces deriving from |
| 137 |
< |
potential energy and constraint forces which are used to guarantee the |
| 138 |
< |
rigidness. However, due to their iterative nature, the SHAKE and |
| 139 |
< |
Rattle algorithms also converge very slowly when the number of |
| 140 |
< |
constraints increases.\cite{Ryckaert1977,Andersen1983} |
| 149 |
> |
Euler angles are a natural choice to describe the rotational degrees |
| 150 |
> |
of freedom. However, due to $1 \over \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 |
> |
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 breakthrough 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 proposed |
| 166 |
< |
to evolve the Hamiltonian system in a constraint manifold by |
| 167 |
< |
iteratively satisfying the orthogonality constraint $Q^T Q = 1$. An |
| 149 |
< |
alternative method using the quaternion representation was developed |
| 150 |
< |
by Omelyan.\cite{Omelyan1998} However, both of these methods are |
| 151 |
< |
iterative and suffer from some related inefficiencies. A symplectic |
| 152 |
< |
Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et |
| 153 |
< |
al.}\cite{Dullweber1997} gets around most of the limitations mentioned |
| 154 |
< |
above and has become the basis for our Langevin integrator. |
| 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 |
+ |
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 |
|
|
| 157 |
– |
\subsection{The Hydrodynamic tensor and Brownian dynamics} |
| 158 |
– |
Combining Brownian dynamics with rigid substructures, one can study |
| 159 |
– |
slow processes in biomolecular systems. Modeling DNA as a chain of |
| 160 |
– |
beads which are subject to harmonic potentials as well as excluded |
| 161 |
– |
volume potentials, Mielke and his coworkers discovered rapid |
| 162 |
– |
superhelical stress generations from the stochastic simulation of twin |
| 163 |
– |
supercoiling DNA with response to induced torques.\cite{Mielke2004} |
| 164 |
– |
Membrane fusion is another key biological process which controls a |
| 165 |
– |
variety of physiological functions, such as release of |
| 166 |
– |
neurotransmitters \textit{etc}. A typical fusion event happens on the |
| 167 |
– |
time scale of a millisecond, which is impractical to study using |
| 168 |
– |
atomistic models with newtonian mechanics. With the help of |
| 169 |
– |
coarse-grained rigid body model and stochastic dynamics, the fusion |
| 170 |
– |
pathways were explored by Noguchi and others.\cite{Noguchi2001,Noguchi2002,Shillcock2005} |
| 171 |
– |
|
| 172 |
– |
Due to the difficulty of numerically integrating anisotropic |
| 173 |
– |
rotational motion, most of the coarse-grained rigid body models are |
| 174 |
– |
treated as spheres, cylinders, ellipsoids or other regular shapes in |
| 175 |
– |
stochastic simulations. In an effort to account for the diffusion |
| 176 |
– |
anisotropy of arbitrarily-shaped particles, Fernandes and Garc\'{i}a |
| 177 |
– |
de la Torre improved the original Brownian dynamics simulation |
| 178 |
– |
algorithm~\cite{Ermak1978,Allison1991} by incorporating a generalized |
| 179 |
– |
$6\times6$ diffusion tensor and introducing a rotational evolution |
| 180 |
– |
scheme consisting of three consecutive rotations.\cite{Fernandes2002} |
| 181 |
– |
Unfortunately, biases are introduced into the system due to the |
| 182 |
– |
arbitrary order of applying the noncommuting rotation |
| 183 |
– |
operators.\cite{Beard2003} Based on the observation the momentum |
| 184 |
– |
relaxation time is much less than the time step, one may ignore the |
| 185 |
– |
inertia in Brownian dynamics. However, the assumption of zero average |
| 186 |
– |
acceleration is not always true for cooperative motion which is common |
| 187 |
– |
in proteins. An inertial Brownian dynamics (IBD) was proposed to |
| 188 |
– |
address this issue by adding an inertial correction |
| 189 |
– |
term.\cite{Beard2000} As a complement to IBD which has a lower bound |
| 190 |
– |
in time step because of the inertial relaxation time, long-time-step |
| 191 |
– |
inertial dynamics (LTID) can be used to investigate the inertial |
| 192 |
– |
behavior of the polymer segments in low friction |
| 193 |
– |
regime.\cite{Beard2000} LTID can also deal with the rotational |
| 194 |
– |
dynamics for nonskew bodies without translation-rotation coupling by |
| 195 |
– |
separating the translation and rotation motion and taking advantage of |
| 196 |
– |
the analytical solution of hydrodynamics properties. However, typical |
| 197 |
– |
nonskew bodies like cylinders and ellipsoids are inadequate to |
| 198 |
– |
represent most complex macromolecular assemblies. |
| 199 |
– |
|
| 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 into |
| 186 |
|
a symplectic rigid body dynamics propagator. In the sections below, |
| 187 |
< |
we review some of the theory of hydrodynamic tensors developed for |
| 188 |
< |
Brownian simulations of rigid multi-particle systems, we then present |
| 189 |
< |
our integration method for a set of generalized Langevin equations of |
| 190 |
< |
motion, and we compare the behavior of the new Langevin integrator to |
| 191 |
< |
dynamical quantities obtained via explicit solvent molecular dynamics. |
| 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 |
|
\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 solute becomes complex. Instead, various |
| 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). |
| 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 |
| 215 |
|
\left( \begin{array}{l} |
| 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} |
| 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). |
| 227 |
|
\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} |
| 228 |
|
For a spherical particle under ``stick'' boundary conditions, the |
| 229 |
|
translational and rotational friction tensors can be calculated from |
| 230 |
< |
Stoke's law, |
| 230 |
> |
Stokes' law, |
| 231 |
|
\begin{equation} |
| 232 |
< |
\Xi^{tt} = \left( {\begin{array}{*{20}c} |
| 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) |
| 236 |
> |
\end{array} \right) |
| 237 |
|
\end{equation} |
| 238 |
|
and |
| 239 |
|
\begin{equation} |
| 240 |
< |
\Xi ^{rr} = \left( {\begin{array}{*{20}c} |
| 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) |
| 244 |
> |
\end{array} \right) |
| 245 |
|
\end{equation} |
| 246 |
|
where $\eta$ is the viscosity of the solvent and $R$ is the |
| 247 |
|
hydrodynamic radius. |
| 249 |
|
Other non-spherical shapes, such as cylinders and ellipsoids, are |
| 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 |
| 270 |
< |
by\cite{Perrin1934,Perrin1936} |
| 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 |
| 255 |
> |
\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1. |
| 256 |
|
\end{equation} |
| 257 |
< |
where the semi-axes are of lengths $a$, $b$, and $c$. Due to the |
| 258 |
< |
complexity of the elliptic integral, only uniaxial ellipsoids, |
| 259 |
< |
{\it i.e.} prolate ($ a \ge b = c$) and oblate ($ a < b = c $), can |
| 260 |
< |
be solved exactly. Introducing an elliptic integral parameter $S$ for |
| 278 |
< |
prolate ellipsoids : |
| 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 ellipsoids: |
| 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 |
< |
one can write down the translational and rotational resistance |
| 269 |
< |
tensors for oblate, |
| 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 |
< |
and prolate, |
| 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 |
< |
ellipsoids. For both spherical and ellipsoidal particles, the |
| 280 |
< |
translation-rotation and rotation-translation coupling tensors are |
| 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}} |
| 302 |
– |
|
| 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 |
| 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. Let us consider a rigid assembly of $N$ beads |
| 300 |
< |
immersed in a continuous medium. Due to hydrodynamic interaction, the |
| 301 |
< |
``net'' velocity of $i$th bead, $v'_i$ is different than its |
| 302 |
< |
unperturbed velocity $v_i$, |
| 303 |
< |
\[ |
| 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 |
|
\[ |
| 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} \\ |
| 456 |
|
\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
| 457 |
|
Consider the Langevin equations of motion in generalized coordinates |
| 458 |
|
\begin{equation} |
| 459 |
< |
M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)} + F_{r,i} (t) |
| 459 |
> |
\mathbf{M}_i \dot \mathbf{V}_i(t) = \mathbf{F}_{s,i}(t) + \mathbf{F}_{f,i}(t) + \mathbf{R}_{i}(t) |
| 460 |
|
\label{LDGeneralizedForm} |
| 461 |
|
\end{equation} |
| 462 |
< |
where $M_i$ is a $6\times6$ generalized diagonal mass (include mass |
| 463 |
< |
and moment of inertial) matrix and $V_i$ is a generalized velocity, |
| 464 |
< |
$V_i = V_i(v_i,\omega _i)$. The right side of |
| 465 |
< |
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in |
| 466 |
< |
lab-fixed frame, systematic force $F_{s,i}$, dissipative force |
| 467 |
< |
$F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the |
| 468 |
< |
system in Newtownian mechanics typically refers to lab-fixed frame, |
| 469 |
< |
it is also convenient to handle the rotation of rigid body in |
| 470 |
< |
body-fixed frame. Thus the friction and random forces are calculated |
| 471 |
< |
in body-fixed frame and converted back to lab-fixed frame by: |
| 472 |
< |
\[ |
| 462 |
> |
where $\mathbf{M}_i$ is a $6\times6$ diagonal mass matrix (which |
| 463 |
> |
includes the rigid body mass and moments of inertia) and $\mathbf{V}_i$ is a |
| 464 |
> |
generalized velocity, $\mathbf{V}_i = |
| 465 |
> |
\left\{\mathbf{v}_i,\mathbf{\omega}_i \right\}$. The right side of |
| 466 |
> |
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a |
| 467 |
> |
system force $\mathbf{F}_{s,i}$, a frictional or dissipative force |
| 468 |
> |
$\mathbf{F}_{f,i}$ and stochastic force $\mathbf{R}_{i}$. While the |
| 469 |
> |
evolution of the system in Newtownian mechanics is typically done in the |
| 470 |
> |
lab-fixed frame, it is convenient to handle the rotation of rigid |
| 471 |
> |
bodies in the body-fixed frame. Thus the friction and random forces are |
| 472 |
> |
calculated in body-fixed frame and converted back to lab-fixed frame |
| 473 |
> |
using the rigid body's rotation matrix ($Q_i$): |
| 474 |
> |
\begin{equation} |
| 475 |
|
\begin{array}{l} |
| 476 |
< |
F_{f,i}^l (t) = Q^T F_{f,i}^b (t), \\ |
| 477 |
< |
F_{r,i}^l (t) = Q^T F_{r,i}^b (t). \\ |
| 476 |
> |
\mathbf{F}_{f,i}(t) = Q_{i}^{T} \mathbf{F}_{f,i}^b (t), \\ |
| 477 |
> |
\mathbf{R}_{i}(t) = Q_{i}^{T} \mathbf{R}_{i}^b (t). \\ |
| 478 |
|
\end{array} |
| 479 |
< |
\] |
| 480 |
< |
Here, the body-fixed friction force $F_{r,i}^b$ is proportional to |
| 481 |
< |
the body-fixed velocity at center of resistance $v_{R,i}^b$ and |
| 482 |
< |
angular velocity $\omega _i$ |
| 479 |
> |
\end{equation} |
| 480 |
> |
Here, the body-fixed friction force $\mathbf{F}_{f,i}^b$ is proportional to |
| 481 |
> |
the body-fixed velocity at the center of resistance $\mathbf{v}_{R,i}^b$ and |
| 482 |
> |
angular velocity $\mathbf{\omega}_i$ |
| 483 |
|
\begin{equation} |
| 484 |
< |
F_{r,i}^b (t) = \left( \begin{array}{l} |
| 485 |
< |
f_{r,i}^b (t) \\ |
| 486 |
< |
\tau _{r,i}^b (t) \\ |
| 487 |
< |
\end{array} \right) = - \left( {\begin{array}{*{20}c} |
| 488 |
< |
{\Xi _{R,t} } & {\Xi _{R,c}^T } \\ |
| 489 |
< |
{\Xi _{R,c} } & {\Xi _{R,r} } \\ |
| 490 |
< |
\end{array}} \right)\left( \begin{array}{l} |
| 491 |
< |
v_{R,i}^b (t) \\ |
| 492 |
< |
\omega _i (t) \\ |
| 484 |
> |
\mathbf{F}_{f,i}^b (t) = \left( \begin{array}{l} |
| 485 |
> |
\mathbf{f}_{f,i}^b (t) \\ |
| 486 |
> |
\mathbf{\tau}_{f,i}^b (t) \\ |
| 487 |
> |
\end{array} \right) = - \left( \begin{array}{*{20}c} |
| 488 |
> |
\Xi_{R,t} & \Xi_{R,c}^T \\ |
| 489 |
> |
\Xi_{R,c} & \Xi_{R,r} \\ |
| 490 |
> |
\end{array} \right)\left( \begin{array}{l} |
| 491 |
> |
\mathbf{v}_{R,i}^b (t) \\ |
| 492 |
> |
\mathbf{\omega}_i (t) \\ |
| 493 |
|
\end{array} \right), |
| 494 |
|
\end{equation} |
| 495 |
< |
while the random force $F_{r,i}^l$ is a Gaussian stochastic variable |
| 495 |
> |
while the random force $\mathbf{R}_{i}^l$ is a Gaussian stochastic variable |
| 496 |
|
with zero mean and variance |
| 497 |
|
\begin{equation} |
| 498 |
< |
\left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle = |
| 499 |
< |
\left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle = |
| 500 |
< |
2k_B T\Xi _R \delta (t - t'). \label{randomForce} |
| 498 |
> |
\left\langle {\mathbf{R}_{i}^l (t) (\mathbf{R}_{i}^l (t'))^T } \right\rangle = |
| 499 |
> |
\left\langle {\mathbf{R}_{i}^b (t) (\mathbf{R}_{i}^b (t'))^T } \right\rangle = |
| 500 |
> |
2 k_B T \Xi_R \delta(t - t'). \label{randomForce} |
| 501 |
|
\end{equation} |
| 502 |
< |
The equation of motion for $v_i$ can be written as |
| 502 |
> |
Once the $6\times6$ resistance tensor at the center of resistance |
| 503 |
> |
($\Xi_R$) is known, obtaining a stochastic vector that has the |
| 504 |
> |
properties in Eq. (\ref{eq:randomForce}) can be done efficiently by |
| 505 |
> |
carrying out a one-time Cholesky decomposition to obtain the square |
| 506 |
> |
root matrix of $\Xi_R$.\cite{SchlickBook} Each time a random force |
| 507 |
> |
vector is needed, a gaussian random vector is generated and then the |
| 508 |
> |
square root matrix is multiplied onto this vector. |
| 509 |
> |
|
| 510 |
> |
The equation of motion for $\mathbf{v}_i$ can be written as |
| 511 |
|
\begin{equation} |
| 512 |
< |
m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) + |
| 513 |
< |
f_{r,i}^l (t) |
| 512 |
> |
m\dot \mathbf{v}_i (t) = \mathbf{f}_{s,i} (t) + \mathbf{f}_{f,i}^l (t) + |
| 513 |
> |
\mathbf{R}_{i}^l (t) |
| 514 |
|
\end{equation} |
| 515 |
|
Since the frictional force is applied at the center of resistance |
| 516 |
|
which generally does not coincide with the center of mass, an extra |
| 517 |
|
torque is exerted at the center of mass. Thus, the net body-fixed |
| 518 |
< |
frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is |
| 518 |
> |
frictional torque at the center of mass, $\tau_{f,i}^b (t)$, is |
| 519 |
|
given by |
| 520 |
|
\begin{equation} |
| 521 |
< |
\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b |
| 521 |
> |
\tau_{f,i}^b \leftarrow \tau_{f,i}^b + \mathbf{r}_{MR} \times \mathbf{f}_{r,i}^b |
| 522 |
|
\end{equation} |
| 523 |
|
where $r_{MR}$ is the vector from the center of mass to the center |
| 524 |
|
of the resistance. Instead of integrating the angular velocity in |
| 525 |
|
lab-fixed frame, we consider the equation of angular momentum in |
| 526 |
|
body-fixed frame |
| 527 |
|
\begin{equation} |
| 528 |
< |
\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t) |
| 537 |
< |
+ \tau _{r,i}^b(t) |
| 528 |
> |
\dot j_i (t) = \tau_{s,i} (t) + \tau_{f,i}^b (t) + \mathbf{R}_{i}^b(t) |
| 529 |
|
\end{equation} |
| 530 |
< |
Embedding the friction terms into force and torque, one can |
| 531 |
< |
integrate the langevin equations of motion for rigid body of |
| 532 |
< |
arbitrary shape in a velocity-Verlet style 2-part algorithm, where |
| 542 |
< |
$h= \delta t$: |
| 530 |
> |
Embedding the friction terms into force and torque, one can integrate |
| 531 |
> |
the Langevin equations of motion for rigid body of arbitrary shape in |
| 532 |
> |
a velocity-Verlet style 2-part algorithm, where $h= \delta t$: |
| 533 |
|
|
| 534 |
|
{\tt moveA:} |
| 535 |
|
\begin{align*} |
