| 39 |
|
\begin{abstract} |
| 40 |
|
We present an algorithm for carrying out Langevin dynamics simulations |
| 41 |
|
on complex rigid bodies by incorporating the hydrodynamic resistance |
| 42 |
< |
tensors for arbitrary shapes into an advanced symplectic integration |
| 42 |
> |
tensors for arbitrary shapes into an advanced rotational integration |
| 43 |
|
scheme. The integrator gives quantitative agreement with both |
| 44 |
|
analytic and approximate hydrodynamic theories for a number of model |
| 45 |
|
rigid bodies, and works well at reproducing the solute dynamical |
| 49 |
|
|
| 50 |
|
\newpage |
| 51 |
|
|
| 52 |
– |
|
| 53 |
– |
|
| 52 |
|
%\narrowtext |
| 53 |
|
|
| 54 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 173 |
|
rigid body) increases.\cite{Ryckaert1977,Andersen1983} |
| 174 |
|
|
| 175 |
|
In order to develop a stable and efficient integration scheme that |
| 176 |
< |
preserves most constants of the motion, symplectic propagators are |
| 177 |
< |
necessary. By introducing a conjugate momentum to the rotation matrix |
| 178 |
< |
${\bf Q}$ and re-formulating Hamilton's equations, a symplectic |
| 179 |
< |
orientational integrator, RSHAKE,\cite{Kol1997} was proposed to evolve |
| 180 |
< |
rigid bodies on a constraint manifold by iteratively satisfying the |
| 181 |
< |
orthogonality constraint ${\bf Q}^T {\bf Q} = 1$. An alternative |
| 182 |
< |
method using the quaternion representation was developed by |
| 183 |
< |
Omelyan.\cite{Omelyan1998} However, both of these methods are |
| 184 |
< |
iterative and suffer from some related inefficiencies. A symplectic |
| 185 |
< |
Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et |
| 186 |
< |
al.}\cite{Dullweber1997} removes most of the limitations mentioned |
| 187 |
< |
above and is therefore the basis for our Langevin integrator. |
| 176 |
> |
preserves most constants of the motion in microcanonical simulations, |
| 177 |
> |
symplectic propagators are necessary. By introducing a conjugate |
| 178 |
> |
momentum to the rotation matrix ${\bf Q}$ and re-formulating |
| 179 |
> |
Hamilton's equations, a symplectic orientational integrator, |
| 180 |
> |
RSHAKE,\cite{Kol1997} was proposed to evolve rigid bodies on a |
| 181 |
> |
constraint manifold by iteratively satisfying the orthogonality |
| 182 |
> |
constraint ${\bf Q}^T {\bf Q} = 1$. An alternative method using the |
| 183 |
> |
quaternion representation was developed by Omelyan.\cite{Omelyan1998} |
| 184 |
> |
However, both of these methods are iterative and suffer from some |
| 185 |
> |
related inefficiencies. A symplectic Lie-Poisson integrator for rigid |
| 186 |
> |
bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes |
| 187 |
> |
most of the limitations mentioned above and is therefore the basis for |
| 188 |
> |
our Langevin integrator. |
| 189 |
|
|
| 190 |
|
The goal of the present work is to develop a Langevin dynamics |
| 191 |
|
algorithm for arbitrary-shaped rigid particles by integrating an |
| 192 |
|
accurate estimate of the friction tensor from hydrodynamics theory |
| 193 |
< |
into a symplectic rigid body dynamics propagator. In the sections |
| 194 |
< |
below, we review some of the theory of hydrodynamic tensors developed |
| 195 |
< |
primarily for Brownian simulations of multi-particle systems, we then |
| 196 |
< |
present our integration method for a set of generalized Langevin |
| 197 |
< |
equations of motion, and we compare the behavior of the new Langevin |
| 198 |
< |
integrator to dynamical quantities obtained via explicit solvent |
| 199 |
< |
molecular dynamics. |
| 193 |
> |
into a stable and efficient rigid body dynamics propagator. In the |
| 194 |
> |
sections below, we review some of the theory of hydrodynamic tensors |
| 195 |
> |
developed primarily for Brownian simulations of multi-particle |
| 196 |
> |
systems, we then present our integration method for a set of |
| 197 |
> |
generalized Langevin equations of motion, and we compare the behavior |
| 198 |
> |
of the new Langevin integrator to dynamical quantities obtained via |
| 199 |
> |
explicit solvent molecular dynamics. |
| 200 |
|
|
| 201 |
|
\subsection{\label{introSection:frictionTensor}The Friction Tensor} |
| 202 |
|
Theoretically, a complete friction kernel for a solute particle can be |
| 231 |
|
\mathbf{\omega} \\ |
| 232 |
|
\end{array} \right). |
| 233 |
|
\end{equation} |
| 234 |
+ |
For an arbitrary body moving in a fluid, Peters has derived a set of |
| 235 |
+ |
fluctuation-dissipation relations for the friction |
| 236 |
+ |
tensors,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk} |
| 237 |
+ |
\begin{eqnarray} |
| 238 |
+ |
\Xi^{tt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf |
| 239 |
+ |
F}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}^2 |
| 240 |
+ |
\right] ds \\ |
| 241 |
+ |
\notag \\ |
| 242 |
+ |
\Xi^{tr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf |
| 243 |
+ |
F}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq} |
| 244 |
+ |
\langle {\bf \tau} \rangle_{eq} \right] ds \\ |
| 245 |
+ |
\notag \\ |
| 246 |
+ |
\Xi^{rt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf |
| 247 |
+ |
\tau}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq} |
| 248 |
+ |
\langle {\bf F} \rangle_{eq} \right] ds \\ |
| 249 |
+ |
\notag \\ |
| 250 |
+ |
\Xi^{rr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf |
| 251 |
+ |
\tau}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}^2 |
| 252 |
+ |
\right] ds |
| 253 |
+ |
\end{eqnarray} |
| 254 |
+ |
In these expressions, the forces (${\bf F}$) and torques (${\bf |
| 255 |
+ |
\tau}$) are those that arise solely from the interactions of the body with |
| 256 |
+ |
the surrounding fluid. For a single solute body in an isotropic fluid, |
| 257 |
+ |
the average forces and torques in these expressions ($\langle {\bf F} |
| 258 |
+ |
\rangle_{eq}$ and $\langle {\bf \tau} \rangle_{eq}$) |
| 259 |
+ |
vanish, and one obtains the simpler force-torque correlation formulae |
| 260 |
+ |
of Nienhuis.\cite{Nienhuis:1970lr} Molecular dynamics simulations with |
| 261 |
+ |
explicit solvent molecules can be used to obtain estimates of the |
| 262 |
+ |
friction tensors with these formulae. In practice, however, one needs |
| 263 |
+ |
relatively long simulations with frequently-stored force and torque |
| 264 |
+ |
information to compute friction tensors, and this becomes |
| 265 |
+ |
prohibitively expensive when there are large numbers of large solute |
| 266 |
+ |
particles. For bodies with simple shapes, there are a number of |
| 267 |
+ |
approximate expressions that allow computation of these tensors |
| 268 |
+ |
without the need for expensive simulations that utilize explicit |
| 269 |
+ |
solvent particles. |
| 270 |
|
|
| 271 |
|
\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} |
| 272 |
< |
For a spherical body under ``stick'' boundary conditions, |
| 273 |
< |
the translational and rotational friction tensors can be calculated |
| 274 |
< |
from Stokes' law, |
| 272 |
> |
For a spherical body under ``stick'' boundary conditions, the |
| 273 |
> |
translational and rotational friction tensors can be estimated from |
| 274 |
> |
Stokes' law, |
| 275 |
|
\begin{equation} |
| 276 |
|
\label{eq:StokesTranslation} |
| 277 |
|
\Xi^{tt} = \left( \begin{array}{*{20}c} |
| 332 |
|
rotation-translation coupling tensors are zero. |
| 333 |
|
|
| 334 |
|
\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}} |
| 335 |
< |
There is no analytical solution for the friction tensor for rigid |
| 336 |
< |
molecules of arbitrary shape. The ellipsoid of revolution and general |
| 337 |
< |
triaxial ellipsoid models have been widely used to approximate the |
| 338 |
< |
hydrodynamic properties of rigid bodies. However, the mapping from all |
| 339 |
< |
possible ellipsoidal spaces ($r$-space) to all possible combinations |
| 340 |
< |
of rotational diffusion coefficients ($D$-space) is not |
| 335 |
> |
Other than the fluctuation dissipation formulae given by |
| 336 |
> |
Peters,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk} there are no |
| 337 |
> |
analytic solutions for the friction tensor for rigid molecules of |
| 338 |
> |
arbitrary shape. The ellipsoid of revolution and general triaxial |
| 339 |
> |
ellipsoid models have been widely used to approximate the hydrodynamic |
| 340 |
> |
properties of rigid bodies. However, the mapping from all possible |
| 341 |
> |
ellipsoidal spaces ($r$-space) to all possible combinations of |
| 342 |
> |
rotational diffusion coefficients ($D$-space) is not |
| 343 |
|
unique.\cite{Wegener1979} Additionally, because there is intrinsic |
| 344 |
|
coupling between translational and rotational motion of {\it skew} |
| 345 |
|
rigid bodies, general ellipsoids are not always suitable for modeling |
| 1196 |
|
|
| 1197 |
|
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
| 1198 |
|
used recently as models for lipid |
| 1199 |
< |
molecules.\cite{SunX._jp0762020,Ayton01} A reference system composed of |
| 1200 |
< |
a single lipid rigid body embedded in a sea of 1929 solvent particles |
| 1201 |
< |
was created and run under a microcanonical ensemble. The resulting |
| 1202 |
< |
viscosity of this mixture was 0.349 centipoise (as estimated using |
| 1203 |
< |
Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties of |
| 1204 |
< |
the lipid rigid body model, we created a rough shell (see |
| 1199 |
> |
molecules.\cite{SunX._jp0762020,Ayton01} A reference system composed |
| 1200 |
> |
of a single lipid rigid body embedded in a sea of 1929 solvent |
| 1201 |
> |
particles was created and run under a microcanonical ensemble. The |
| 1202 |
> |
resulting viscosity of this mixture was 0.349 centipoise (as estimated |
| 1203 |
> |
using Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties |
| 1204 |
> |
of the lipid rigid body model, we created a rough shell (see |
| 1205 |
|
Fig.~\ref{fig:roughShell}), in which the lipid is represented as a |
| 1206 |
|
``shell'' made of 3550 identical beads (0.25 \AA\ in diameter) |
| 1207 |
< |
distributed on the surface. Applying the procedure described in |
| 1208 |
< |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
| 1207 |
> |
distributed on the surface. Applying the procedure described by |
| 1208 |
> |
Eq. (\ref{introEquation:ResistanceTensorArbitraryOrigin}), we |
| 1209 |
|
identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46 |
| 1210 |
|
\AA). |
| 1211 |
|
|
| 1346 |
|
|
| 1347 |
|
We have presented a new algorithm for carrying out Langevin dynamics |
| 1348 |
|
simulations on complex rigid bodies by incorporating the hydrodynamic |
| 1349 |
< |
resistance tensors for arbitrary shapes into an advanced symplectic |
| 1349 |
> |
resistance tensors for arbitrary shapes into a stable and efficient |
| 1350 |
|
integration scheme. The integrator gives quantitative agreement with |
| 1351 |
|
both analytic and approximate hydrodynamic theories, and works |
| 1352 |
|
reasonably well at reproducing the solute dynamical properties |
