| 1 |
|
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
| 2 |
< |
%\documentclass[aps,prb,preprint]{revtex4} |
| 2 |
> |
%\documentclass[aps,prb,preprint,amsmath,amssymb,endfloats]{revtex4} |
| 3 |
|
\documentclass[11pt]{article} |
| 4 |
– |
\usepackage{endfloat} |
| 4 |
|
\usepackage{amsmath} |
| 5 |
|
\usepackage{amssymb} |
| 6 |
|
\usepackage{times} |
| 7 |
|
\usepackage{mathptmx} |
| 8 |
|
\usepackage{setspace} |
| 9 |
< |
\usepackage{tabularx} |
| 9 |
> |
\usepackage{endfloat} |
| 10 |
> |
%\usepackage{tabularx} |
| 11 |
|
\usepackage{graphicx} |
| 12 |
< |
\usepackage{booktabs} |
| 13 |
< |
\usepackage{bibentry} |
| 14 |
< |
\usepackage{mathrsfs} |
| 12 |
> |
%\usepackage{booktabs} |
| 13 |
> |
%\usepackage{bibentry} |
| 14 |
> |
%\usepackage{mathrsfs} |
| 15 |
|
\usepackage[ref]{overcite} |
| 16 |
|
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm |
| 17 |
|
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight |
| 18 |
|
9.0in \textwidth 6.5in \brokenpenalty=10000 |
| 19 |
< |
\renewcommand{\baselinestretch}{1.2} |
| 19 |
> |
|
| 20 |
|
\renewcommand\citemid{\ } % no comma in optional referenc note |
| 21 |
|
|
| 22 |
|
\begin{document} |
| 23 |
|
|
| 24 |
< |
\title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape } |
| 24 |
> |
\title{Langevin dynamics for rigid bodies of arbitrary shape} |
| 25 |
|
|
| 26 |
< |
\author{Xiuquan Sun, Teng Lin and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 27 |
< |
gezelter@nd.edu} \\ |
| 28 |
< |
Department of Chemistry and Biochemistry\\ |
| 26 |
> |
\author{Xiuquan Sun, Teng Lin and J. Daniel |
| 27 |
> |
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ |
| 28 |
> |
Department of Chemistry and Biochemistry,\\ |
| 29 |
|
University of Notre Dame\\ |
| 30 |
|
Notre Dame, Indiana 46556} |
| 31 |
|
|
| 32 |
|
\date{\today} |
| 33 |
|
|
| 34 |
– |
\maketitle \doublespacing |
| 34 |
|
|
| 35 |
< |
\begin{abstract} |
| 35 |
> |
\maketitle |
| 36 |
|
|
| 37 |
+ |
|
| 38 |
+ |
|
| 39 |
+ |
\begin{abstract} |
| 40 |
+ |
We present an algorithm for carrying out Langevin dynamics simulations |
| 41 |
+ |
on complex rigid bodies by incorporating the hydrodynamic resistance |
| 42 |
+ |
tensors for arbitrary shapes into an advanced symplectic integration |
| 43 |
+ |
scheme. The integrator gives quantitative agreement with both |
| 44 |
+ |
analytic and approximate hydrodynamic theories for a number of model |
| 45 |
+ |
rigid bodies, and works well at reproducing the solute dynamical |
| 46 |
+ |
properties (diffusion constants, and orientational relaxation times) |
| 47 |
+ |
obtained from explicitly-solvated simulations. |
| 48 |
|
\end{abstract} |
| 49 |
|
|
| 50 |
|
\newpage |
| 51 |
|
|
| 52 |
+ |
|
| 53 |
+ |
|
| 54 |
|
%\narrowtext |
| 55 |
|
|
| 56 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 57 |
|
% BODY OF TEXT |
| 58 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 59 |
|
|
| 60 |
+ |
\begin{doublespace} |
| 61 |
+ |
|
| 62 |
|
\section{Introduction} |
| 63 |
|
|
| 64 |
|
%applications of langevin dynamics |
| 65 |
< |
Langevin dynamics, which mimics a simple heat bath with stochastic and |
| 65 |
> |
Langevin dynamics, which mimics a heat bath using both stochastic and |
| 66 |
|
dissipative forces, has been applied in a variety of situations as an |
| 67 |
|
alternative to molecular dynamics with explicit solvent molecules. |
| 68 |
|
The stochastic treatment of the solvent allows the use of simulations |
| 69 |
< |
with substantially longer time and length scales. In general, the |
| 69 |
> |
with substantially longer time and length scales. In general, the |
| 70 |
|
dynamic and structural properties obtained from Langevin simulations |
| 71 |
|
agree quite well with similar properties obtained from explicit |
| 72 |
|
solvent simulations. |
| 75 |
|
study of met-enkephalin in which Langevin simulations predicted |
| 76 |
|
dynamical properties that were largely in agreement with explicit |
| 77 |
|
solvent simulations.\cite{Shen2002} By applying Langevin dynamics with |
| 78 |
< |
the UNRES model, Liow and his coworkers suggest that protein folding |
| 78 |
> |
the UNRES model, Liwo and his coworkers suggest that protein folding |
| 79 |
|
pathways can be explored within a reasonable amount of |
| 80 |
|
time.\cite{Liwo2005} |
| 81 |
|
|
| 90 |
|
folding/unfolding studies and discovered a higher free energy barrier |
| 91 |
|
between the native and denatured states.\cite{HuseyinKaya07012005} |
| 92 |
|
|
| 93 |
< |
Because of its stability against noise, Langevin dynamics has also |
| 80 |
< |
proven useful for studying remagnetization processes in various |
| 81 |
< |
systems.\cite{Palacios1998,Berkov2002,Denisov2003} [Check: For |
| 82 |
< |
instance, the oscillation power spectrum of nanoparticles from |
| 83 |
< |
Langevin dynamics has the same peak frequencies for different wave |
| 84 |
< |
vectors, which recovers the property of magnetic excitations in small |
| 85 |
< |
finite structures.\cite{Berkov2005a}] |
| 86 |
< |
|
| 87 |
< |
In typical LD simulations, the friction and random forces on |
| 93 |
> |
In typical LD simulations, the friction and random ($f_r$) forces on |
| 94 |
|
individual atoms are taken from Stokes' law, |
| 95 |
|
\begin{eqnarray} |
| 96 |
< |
m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\ |
| 97 |
< |
\langle R(t) \rangle & = & 0 \\ |
| 98 |
< |
\langle R(t) R(t') \rangle & = & 2 k_B T \xi m \delta(t - t') |
| 96 |
> |
m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + f_r(t) \notag \\ |
| 97 |
> |
\langle f_r(t) \rangle & = & 0 \\ |
| 98 |
> |
\langle f_r(t) f_r(t') \rangle & = & 2 k_B T \xi m \delta(t - t') \notag |
| 99 |
|
\end{eqnarray} |
| 100 |
< |
where $\xi \approx 6 \pi \eta a$. Here $\eta$ is the viscosity of the |
| 101 |
< |
implicit solvent, and $a$ is the hydrodynamic radius of the atom. |
| 100 |
> |
where $\xi \approx 6 \pi \eta \rho$. Here $\eta$ is the viscosity of the |
| 101 |
> |
implicit solvent, and $\rho$ is the hydrodynamic radius of the atom. |
| 102 |
|
|
| 103 |
|
The use of rigid substructures,\cite{Chun:2000fj} |
| 104 |
< |
coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunGezelter08} |
| 105 |
< |
and ellipsoidal representations of protein side chains~\cite{Fogolari:1996lr} |
| 106 |
< |
has made the use of the Stokes-Einstein approximation problematic. A |
| 107 |
< |
rigid substructure moves as a single unit with orientational as well |
| 108 |
< |
as translational degrees of freedom. This requires a more general |
| 109 |
< |
treatment of the hydrodynamics than the spherical approximation |
| 110 |
< |
provides. The atoms involved in a rigid or coarse-grained structure |
| 111 |
< |
should properly have solvent-mediated interactions with each |
| 112 |
< |
other. The theory of interactions {\it between} bodies moving through |
| 113 |
< |
a fluid has been developed over the past century and has been applied |
| 114 |
< |
to simulations of Brownian |
| 115 |
< |
motion.\cite{FIXMAN:1986lr,Ramachandran1996} |
| 104 |
> |
coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunX._jp0762020} |
| 105 |
> |
and ellipsoidal representations of protein side |
| 106 |
> |
chains~\cite{Fogolari:1996lr} has made the use of the Stokes-Einstein |
| 107 |
> |
approximation problematic. A rigid substructure moves as a single |
| 108 |
> |
unit with orientational as well as translational degrees of freedom. |
| 109 |
> |
This requires a more general treatment of the hydrodynamics than the |
| 110 |
> |
spherical approximation provides. Also, the atoms involved in a rigid |
| 111 |
> |
or coarse-grained structure have solvent-mediated interactions with |
| 112 |
> |
each other, and these interactions are ignored if all atoms are |
| 113 |
> |
treated as separate spherical particles. The theory of interactions |
| 114 |
> |
{\it between} bodies moving through a fluid has been developed over |
| 115 |
> |
the past century and has been applied to simulations of Brownian |
| 116 |
> |
motion.\cite{FIXMAN:1986lr,Ramachandran1996} |
| 117 |
|
|
| 118 |
< |
In order to account for the diffusion anisotropy of arbitrarily-shaped |
| 119 |
< |
particles, Fernandes and Garc\'{i}a de la Torre improved the original |
| 120 |
< |
Brownian dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by |
| 118 |
> |
In order to account for the diffusion anisotropy of complex shapes, |
| 119 |
> |
Fernandes and Garc\'{i}a de la Torre improved an earlier Brownian |
| 120 |
> |
dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by |
| 121 |
|
incorporating a generalized $6\times6$ diffusion tensor and |
| 122 |
|
introducing a rotational evolution scheme consisting of three |
| 123 |
|
consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are |
| 128 |
|
assumption of zero average acceleration is not always true for |
| 129 |
|
cooperative motion which is common in proteins. An inertial Brownian |
| 130 |
|
dynamics (IBD) was proposed to address this issue by adding an |
| 131 |
< |
inertial correction term.\cite{Beard2000} As a complement to IBD which |
| 132 |
< |
has a lower bound in time step because of the inertial relaxation |
| 133 |
< |
time, long-time-step inertial dynamics (LTID) can be used to |
| 134 |
< |
investigate the inertial behavior of linked polymer segments in a low |
| 135 |
< |
friction regime.\cite{Beard2000} LTID can also deal with the |
| 131 |
> |
inertial correction term.\cite{Beard2000} As a complement to IBD, |
| 132 |
> |
which has a lower bound in time step because of the inertial |
| 133 |
> |
relaxation time, long-time-step inertial dynamics (LTID) can be used |
| 134 |
> |
to investigate the inertial behavior of linked polymer segments in a |
| 135 |
> |
low friction regime.\cite{Beard2000} LTID can also deal with the |
| 136 |
|
rotational dynamics for nonskew bodies without translation-rotation |
| 137 |
|
coupling by separating the translation and rotation motion and taking |
| 138 |
< |
advantage of the analytical solution of hydrodynamics |
| 138 |
> |
advantage of the analytical solution of hydrodynamic |
| 139 |
|
properties. However, typical nonskew bodies like cylinders and |
| 140 |
|
ellipsoids are inadequate to represent most complex macromolecular |
| 141 |
< |
assemblies. There is therefore a need for incorporating the |
| 142 |
< |
hydrodynamics of complex (and potentially skew) rigid bodies in the |
| 143 |
< |
library of methods available for performing Langevin simulations. |
| 141 |
> |
assemblies. Therefore, the goal of this work is to adapt some of the |
| 142 |
> |
hydrodynamic methodologies developed to treat Brownian motion of |
| 143 |
> |
complex assemblies into a Langevin integrator for rigid bodies with |
| 144 |
> |
arbitrary shapes. |
| 145 |
|
|
| 146 |
|
\subsection{Rigid Body Dynamics} |
| 147 |
|
Rigid bodies are frequently involved in the modeling of large |
| 148 |
|
collections of particles that move as a single unit. In molecular |
| 149 |
|
simulations, rigid bodies have been used to simplify protein-protein |
| 150 |
|
docking,\cite{Gray2003} and lipid bilayer |
| 151 |
< |
simulations.\cite{SunGezelter08} Many of the water models in common |
| 151 |
> |
simulations.\cite{SunX._jp0762020} Many of the water models in common |
| 152 |
|
use are also rigid-body |
| 153 |
|
models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are |
| 154 |
< |
typically evolved using constraints rather than rigid body equations |
| 155 |
< |
of motion. |
| 154 |
> |
typically evolved in molecular dynamics simulations using constraints |
| 155 |
> |
rather than rigid body equations of motion. |
| 156 |
|
|
| 157 |
|
Euler angles are a natural choice to describe the rotational degrees |
| 158 |
|
of freedom. However, due to $\frac{1}{\sin \theta}$ singularities, the |
| 177 |
|
In order to develop a stable and efficient integration scheme that |
| 178 |
|
preserves most constants of the motion, symplectic propagators are |
| 179 |
|
necessary. By introducing a conjugate momentum to the rotation matrix |
| 180 |
< |
$Q$ and re-formulating Hamilton's equations, a symplectic |
| 180 |
> |
${\bf Q}$ and re-formulating Hamilton's equations, a symplectic |
| 181 |
|
orientational integrator, RSHAKE,\cite{Kol1997} was proposed to evolve |
| 182 |
|
rigid bodies on a constraint manifold by iteratively satisfying the |
| 183 |
< |
orthogonality constraint $Q^T Q = 1$. An alternative method using the |
| 184 |
< |
quaternion representation was developed by Omelyan.\cite{Omelyan1998} |
| 185 |
< |
However, both of these methods are iterative and suffer from some |
| 186 |
< |
related inefficiencies. A symplectic Lie-Poisson integrator for rigid |
| 187 |
< |
bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes |
| 188 |
< |
most of the limitations mentioned above and is therefore the basis for |
| 189 |
< |
our Langevin integrator. |
| 183 |
> |
orthogonality constraint ${\bf Q}^T {\bf Q} = 1$. An alternative |
| 184 |
> |
method using the quaternion representation was developed by |
| 185 |
> |
Omelyan.\cite{Omelyan1998} However, both of these methods are |
| 186 |
> |
iterative and suffer from some related inefficiencies. A symplectic |
| 187 |
> |
Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et |
| 188 |
> |
al.}\cite{Dullweber1997} removes most of the limitations mentioned |
| 189 |
> |
above and is therefore the basis for our Langevin integrator. |
| 190 |
|
|
| 191 |
|
The goal of the present work is to develop a Langevin dynamics |
| 192 |
|
algorithm for arbitrary-shaped rigid particles by integrating an |
| 234 |
|
\end{equation} |
| 235 |
|
|
| 236 |
|
\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} |
| 237 |
< |
For a spherical body under ``stick'' boundary conditions, the |
| 238 |
< |
translational and rotational friction tensors can be calculated from |
| 239 |
< |
Stokes' law,\cite{stokes} |
| 237 |
> |
For a spherical body under ``stick'' boundary conditions, |
| 238 |
> |
the translational and rotational friction tensors can be calculated |
| 239 |
> |
from Stokes' law, |
| 240 |
|
\begin{equation} |
| 241 |
+ |
\label{eq:StokesTranslation} |
| 242 |
|
\Xi^{tt} = \left( \begin{array}{*{20}c} |
| 243 |
< |
{6\pi \eta R} & 0 & 0 \\ |
| 244 |
< |
0 & {6\pi \eta R} & 0 \\ |
| 245 |
< |
0 & 0 & {6\pi \eta R} \\ |
| 243 |
> |
{6\pi \eta \rho} & 0 & 0 \\ |
| 244 |
> |
0 & {6\pi \eta \rho} & 0 \\ |
| 245 |
> |
0 & 0 & {6\pi \eta \rho} \\ |
| 246 |
|
\end{array} \right) |
| 247 |
|
\end{equation} |
| 248 |
|
and |
| 249 |
|
\begin{equation} |
| 250 |
+ |
\label{eq:StokesRotation} |
| 251 |
|
\Xi^{rr} = \left( \begin{array}{*{20}c} |
| 252 |
< |
{8\pi \eta R^3 } & 0 & 0 \\ |
| 253 |
< |
0 & {8\pi \eta R^3 } & 0 \\ |
| 254 |
< |
0 & 0 & {8\pi \eta R^3 } \\ |
| 252 |
> |
{8\pi \eta \rho^3 } & 0 & 0 \\ |
| 253 |
> |
0 & {8\pi \eta \rho^3 } & 0 \\ |
| 254 |
> |
0 & 0 & {8\pi \eta \rho^3 } \\ |
| 255 |
|
\end{array} \right) |
| 256 |
|
\end{equation} |
| 257 |
< |
where $\eta$ is the viscosity of the solvent and $R$ is the |
| 258 |
< |
hydrodynamic radius. |
| 257 |
> |
where $\eta$ is the viscosity of the solvent and $\rho$ is the |
| 258 |
> |
hydrodynamic radius. The presence of the rotational resistance tensor |
| 259 |
> |
implies that the spherical body has internal structure and |
| 260 |
> |
orientational degrees of freedom that must be propagated in time. For |
| 261 |
> |
non-structured spherical bodies (i.e. the atoms in a traditional |
| 262 |
> |
molecular dynamics simulation) these degrees of freedom do not exist. |
| 263 |
|
|
| 264 |
|
Other non-spherical shapes, such as cylinders and ellipsoids, are |
| 265 |
|
widely used as references for developing new hydrodynamic theories, |
| 330 |
|
The frictional force felt by the $i^\mathrm{th}$ bead is proportional to |
| 331 |
|
its net velocity |
| 332 |
|
\begin{equation} |
| 333 |
< |
{\bf F}_i = \zeta_i {\bf v}_i - \zeta _i \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }. |
| 333 |
> |
{\bf F}_i = \xi_i {\bf v}_i - \xi_i \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }. |
| 334 |
|
\label{introEquation:tensorExpression} |
| 335 |
|
\end{equation} |
| 336 |
|
Eq. (\ref{introEquation:tensorExpression}) defines the two-point |
| 338 |
|
solutions to this equation, including the simple solution given by |
| 339 |
|
Oseen and Burgers in 1930 for two beads of identical radius. A second |
| 340 |
|
order expression for beads of different hydrodynamic radii was |
| 341 |
< |
introduced by Rotne and Prager\cite{Rotne1969} and improved by |
| 341 |
> |
introduced by Rotne and Prager,\cite{Rotne1969} and improved by |
| 342 |
|
Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977} |
| 343 |
|
\begin{equation} |
| 344 |
|
{\bf T}_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {{\bf I} + |
| 345 |
< |
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\sigma |
| 346 |
< |
_i^2 + \sigma _j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} - |
| 345 |
> |
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\rho |
| 346 |
> |
_i^2 + \rho_j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} - |
| 347 |
|
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. |
| 348 |
|
\label{introEquation:RPTensorNonOverlapped} |
| 349 |
|
\end{equation} |
| 350 |
|
Here ${\bf R}_{ij}$ is the distance vector between beads $i$ and $j$. Both |
| 351 |
|
the Oseen-Burgers tensor and |
| 352 |
|
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption that |
| 353 |
< |
the beads do not overlap ($R_{ij} \ge \sigma _i + \sigma _j$). |
| 353 |
> |
the beads do not overlap ($R_{ij} \ge \rho_i + \rho_j$). |
| 354 |
|
|
| 355 |
|
To calculate the resistance tensor for a body represented as the union |
| 356 |
|
of many non-overlapping beads, we first pick an arbitrary origin $O$ |
| 403 |
|
additive correction uses the solvent viscosity ($\eta$) as well as the |
| 404 |
|
total volume of the beads that contribute to the hydrodynamic model, |
| 405 |
|
\begin{equation} |
| 406 |
< |
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3, |
| 406 |
> |
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3, |
| 407 |
|
\end{equation} |
| 408 |
< |
where $\sigma_i$ is the radius of bead $i$. This correction term was |
| 408 |
> |
where $\rho_i$ is the radius of bead $i$. This correction term was |
| 409 |
|
rigorously tested and compared with the analytical results for |
| 410 |
|
two-sphere and ellipsoidal systems by Garc\'{i}a de la Torre and |
| 411 |
|
Rodes.\cite{Torre:1983lr} |
| 468 |
|
to be the bead diameter, so that adjacent beads are touching, but do |
| 469 |
|
not overlap. To make a shape corresponding to the rigid structure, |
| 470 |
|
beads that sit on lattice sites that are outside the van der Waals |
| 471 |
< |
radii of any atoms comprising the rigid body are excluded from the |
| 472 |
< |
calculation. |
| 471 |
> |
radii of all of the atoms comprising the rigid body are excluded from |
| 472 |
> |
the calculation. |
| 473 |
|
|
| 474 |
|
For large structures, most of the beads will be deep within the rigid |
| 475 |
|
body and will not contribute to the hydrodynamic tensor. In the {\it |
| 481 |
|
truncation can still produce bead assemblies with thousands of |
| 482 |
|
members. |
| 483 |
|
|
| 484 |
< |
If all of the atoms comprising the rigid substructure are spherical |
| 485 |
< |
and non-overlapping, the tensor in |
| 484 |
> |
If all of the {\it atoms} comprising the rigid substructure are |
| 485 |
> |
spherical and non-overlapping, the tensor in |
| 486 |
|
Eq.~(\ref{introEquation:RPTensorNonOverlapped}) may be used directly |
| 487 |
|
using the atoms themselves as the hydrodynamic beads. This is a |
| 488 |
|
variant of the {\it bead model} approach of Carrasco and Garc\'{i}a de |
| 489 |
< |
la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$ matrix |
| 490 |
< |
can be quite small, and the calculation of the hydrodynamic tensor is |
| 491 |
< |
straightforward. |
| 489 |
> |
la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$ |
| 490 |
> |
matrix can be quite small, and the calculation of the hydrodynamic |
| 491 |
> |
tensor is straightforward. |
| 492 |
|
|
| 493 |
|
In general, the inversion of the ${\bf B}$ matrix is the most |
| 494 |
|
computationally demanding task. This inversion is done only once for |
| 495 |
< |
each type of rigid structure. We have been using straightforward |
| 496 |
< |
LU-decomposition to solve the linear system and obtain the elements of |
| 497 |
< |
${\bf C}$. Once ${\bf C}$ has been obtained, the location of the |
| 495 |
> |
each type of rigid structure. We have used straightforward |
| 496 |
> |
LU-decomposition to solve the linear system and to obtain the elements |
| 497 |
> |
of ${\bf C}$. Once ${\bf C}$ has been obtained, the location of the |
| 498 |
|
center of resistance ($R$) is found and the resistance tensor at this |
| 499 |
|
point is calculated. The $3 \times 1$ vector giving the location of |
| 500 |
|
the rigid body's center of resistance and the $6 \times 6$ resistance |
| 501 |
< |
tensor are stored for use in the Langevin dynamics calculation. Note |
| 502 |
< |
that these quantities depend on solvent viscosity and temperature and |
| 503 |
< |
must be recomputed if different simulation conditions are required. |
| 501 |
> |
tensor are then stored for use in the Langevin dynamics calculation. |
| 502 |
> |
These quantities depend on solvent viscosity and temperature and must |
| 503 |
> |
be recomputed if different simulation conditions are required. |
| 504 |
|
|
| 505 |
|
\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
| 506 |
|
|
| 516 |
|
${\bf V} = |
| 517 |
|
\left\{{\bf v},{\bf \omega}\right\}$. The right side of |
| 518 |
|
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a |
| 519 |
< |
system force (${\bf F}_{s}$), a frictional or dissipative force |
| 520 |
< |
(${\bf F}_{f}$) and stochastic force (${\bf F}_{r}$). While the |
| 521 |
< |
evolution of the system in Newtonian mechanics is typically done in |
| 522 |
< |
the lab frame, it is convenient to handle the dynamics of rigid bodies |
| 523 |
< |
in body-fixed frames. Thus the friction and random forces on each |
| 519 |
> |
system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf |
| 520 |
> |
F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution |
| 521 |
> |
of the system in Newtonian mechanics is typically done in the lab |
| 522 |
> |
frame, it is convenient to handle the dynamics of rigid bodies in |
| 523 |
> |
body-fixed frames. Thus the friction and random forces on each |
| 524 |
|
substructure are calculated in a body-fixed frame and may converted |
| 525 |
|
back to the lab frame using that substructure's rotation matrix (${\bf |
| 526 |
|
Q}$): |
| 568 |
|
\Xi_R = {\bf S} {\bf S}^{T}, |
| 569 |
|
\label{eq:Cholesky} |
| 570 |
|
\end{equation} |
| 571 |
< |
where ${\bf S}$ is a lower triangular matrix.\cite{SchlickBook} A |
| 571 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
| 572 |
|
vector with the statistics required for the random force can then be |
| 573 |
|
obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which |
| 574 |
|
has elements chosen from a Gaussian distribution, such that: |
| 599 |
|
frame, we consider the equation of motion for the angular momentum |
| 600 |
|
(${\bf j}$) in the body-fixed frame |
| 601 |
|
\begin{equation} |
| 602 |
< |
\dot{\bf j}(t) = \tau^{~b}(t) |
| 602 |
> |
\frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t) |
| 603 |
|
\end{equation} |
| 604 |
|
Embedding the friction and random forces into the the total force and |
| 605 |
|
torque, one can integrate the Langevin equations of motion for a rigid |
| 769 |
|
We performed reference microcanonical simulations with explicit |
| 770 |
|
solvents for each of the different model system. In each case there |
| 771 |
|
was one solute model and 1929 solvent molecules present in the |
| 772 |
< |
simulation box. All simulations were equilibrated using a |
| 772 |
> |
simulation box. All simulations were equilibrated for 5 ns using a |
| 773 |
|
constant-pressure and temperature integrator with target values of 300 |
| 774 |
|
K for the temperature and 1 atm for pressure. Following this stage, |
| 775 |
< |
further equilibration and sampling was done in a microcanonical |
| 776 |
< |
ensemble. Since the model bodies are typically quite massive, we were |
| 777 |
< |
able to use a time step of 25 fs. |
| 775 |
> |
further equilibration (5 ns) and sampling (10 ns) was done in a |
| 776 |
> |
microcanonical ensemble. Since the model bodies are typically quite |
| 777 |
> |
massive, we were able to use a time step of 25 fs. |
| 778 |
|
|
| 779 |
|
The model systems studied used both Lennard-Jones spheres as well as |
| 780 |
|
uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne |
| 781 |
|
potential was a single site model for the interactions of rigid |
| 782 |
< |
ellipsoidal molecules.\cite{Gay81} It can be thought of as a |
| 782 |
> |
ellipsoidal molecules.\cite{Gay1981} It can be thought of as a |
| 783 |
|
modification of the Gaussian overlap model originally described by |
| 784 |
|
Berne and Pechukas.\cite{Berne72} The potential is constructed in the |
| 785 |
|
familiar form of the Lennard-Jones function using |
| 808 |
|
Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e / |
| 809 |
|
\epsilon^s$, describes the ratio between the well depths in the {\it |
| 810 |
|
end-to-end} and side-by-side configurations. Details of the potential |
| 811 |
< |
are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an |
| 811 |
> |
are given elsewhere,\cite{Luckhurst90,Golubkov06,SunX._jp0762020} and an |
| 812 |
|
excellent overview of the computational methods that can be used to |
| 813 |
|
efficiently compute forces and torques for this potential can be found |
| 814 |
|
in Ref. \citen{Golubkov06} |
| 843 |
|
\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}. |
| 844 |
|
\label{eq:shear} |
| 845 |
|
\end{equation} |
| 846 |
< |
A similar form exists for the bulk viscosity |
| 847 |
< |
\begin{equation} |
| 834 |
< |
\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left( |
| 835 |
< |
\int_{t_0}^{t_0 + t} |
| 836 |
< |
\left(P\left(t'\right)-\left\langle P \right\rangle \right)dt' |
| 837 |
< |
\right)^2 \right\rangle_{t_0}. |
| 838 |
< |
\end{equation} |
| 839 |
< |
Alternatively, the shear viscosity can also be calculated using a |
| 840 |
< |
Green-Kubo formula with the off-diagonal pressure tensor correlation function, |
| 841 |
< |
\begin{equation} |
| 842 |
< |
\eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0 |
| 843 |
< |
+ t) \right\rangle_{t_0} dt, |
| 844 |
< |
\end{equation} |
| 845 |
< |
although this method converges extremely slowly and is not practical |
| 846 |
< |
for obtaining viscosities from molecular dynamics simulations. |
| 846 |
> |
which converges much more rapidly in molecular dynamics simulations |
| 847 |
> |
than the traditional Green-Kubo formula. |
| 848 |
|
|
| 849 |
|
The Langevin dynamics for the different model systems were performed |
| 850 |
|
at the same temperature as the average temperature of the |
| 870 |
|
compute the diffusive behavior for motion parallel to each body-fixed |
| 871 |
|
axis by projecting the displacement of the particle onto the |
| 872 |
|
body-fixed reference frame at $t=0$. With an isotropic solvent, as we |
| 873 |
< |
have used in this study, there are differences between the three |
| 874 |
< |
diffusion constants, but these must converge to the same value at |
| 875 |
< |
longer times. Translational diffusion constants for the different |
| 876 |
< |
shaped models are shown in table \ref{tab:translation}. |
| 873 |
> |
have used in this study, there may be differences between the three |
| 874 |
> |
diffusion constants at short times, but these must converge to the |
| 875 |
> |
same value at longer times. Translational diffusion constants for the |
| 876 |
> |
different shaped models are shown in table \ref{tab:translation}. |
| 877 |
|
|
| 878 |
|
In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational |
| 879 |
|
diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object |
| 953 |
|
an arbitrary value of 0.8 kcal/mol. |
| 954 |
|
|
| 955 |
|
The Stokes-Einstein behavior of large spherical particles in |
| 956 |
< |
hydrodynamic flows is well known, giving translational friction |
| 957 |
< |
coefficients of $6 \pi \eta R$ (stick boundary conditions) and |
| 958 |
< |
rotational friction coefficients of $8 \pi \eta R^3$. Recently, |
| 959 |
< |
Schmidt and Skinner have computed the behavior of spherical tag |
| 960 |
< |
particles in molecular dynamics simulations, and have shown that {\it |
| 961 |
< |
slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more |
| 962 |
< |
appropriate for molecule-sized spheres embedded in a sea of spherical |
| 963 |
< |
solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
| 956 |
> |
hydrodynamic flows with ``stick'' boundary conditions is well known, |
| 957 |
> |
and is given in Eqs. (\ref{eq:StokesTranslation}) and |
| 958 |
> |
(\ref{eq:StokesRotation}). Recently, Schmidt and Skinner have |
| 959 |
> |
computed the behavior of spherical tag particles in molecular dynamics |
| 960 |
> |
simulations, and have shown that {\it slip} boundary conditions |
| 961 |
> |
($\Xi_{tt} = 4 \pi \eta \rho$) may be more appropriate for |
| 962 |
> |
molecule-sized spheres embedded in a sea of spherical solvent |
| 963 |
> |
particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
| 964 |
|
|
| 965 |
|
Our simulation results show similar behavior to the behavior observed |
| 966 |
|
by Schmidt and Skinner. The diffusion constant obtained from our |
| 985 |
|
can be combined to give a single translational diffusion |
| 986 |
|
constant,\cite{Berne90} |
| 987 |
|
\begin{equation} |
| 988 |
< |
D = \frac{k_B T}{6 \pi \eta a} G(\rho), |
| 988 |
> |
D = \frac{k_B T}{6 \pi \eta a} G(s), |
| 989 |
|
\label{Dperrin} |
| 990 |
|
\end{equation} |
| 991 |
|
as well as a single rotational diffusion coefficient, |
| 992 |
|
\begin{equation} |
| 993 |
< |
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2) |
| 994 |
< |
G(\rho) - 1}{1 - \rho^4} \right\}. |
| 993 |
> |
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - s^2) |
| 994 |
> |
G(s) - 1}{1 - s^4} \right\}. |
| 995 |
|
\label{ThetaPerrin} |
| 996 |
|
\end{equation} |
| 997 |
< |
In these expressions, $G(\rho)$ is a function of the axial ratio |
| 998 |
< |
($\rho = b / a$), which for prolate ellipsoids, is |
| 997 |
> |
In these expressions, $G(s)$ is a function of the axial ratio |
| 998 |
> |
($s = b / a$), which for prolate ellipsoids, is |
| 999 |
|
\begin{equation} |
| 1000 |
< |
G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 - |
| 1000 |
< |
\rho^2)^{1/2}}{\rho} \right\} |
| 1000 |
> |
G(s) = (1- s^2)^{-1/2} \ln \left\{ \frac{1 + (1 - s^2)^{1/2}}{s} \right\} |
| 1001 |
|
\label{GPerrin} |
| 1002 |
|
\end{equation} |
| 1003 |
|
Again, there is some uncertainty about the correct boundary conditions |
| 1022 |
|
exact treatment of the diffusion tensor as well as the rough-shell |
| 1023 |
|
model for the ellipsoid. |
| 1024 |
|
|
| 1025 |
< |
The translational diffusion constants from the microcanonical simulations |
| 1026 |
< |
agree well with the predictions of the Perrin model, although the rotational |
| 1027 |
< |
correlation times are a factor of 2 shorter than expected from hydrodynamic |
| 1028 |
< |
theory. One explanation for the slower rotation |
| 1029 |
< |
of explicitly-solvated ellipsoids is the possibility that solute-solvent |
| 1030 |
< |
collisions happen at both ends of the solute whenever the principal |
| 1031 |
< |
axis of the ellipsoid is turning. In the upper portion of figure |
| 1032 |
< |
\ref{fig:explanation} we sketch a physical picture of this explanation. |
| 1033 |
< |
Since our Langevin integrator is providing nearly quantitative agreement with |
| 1034 |
< |
the Perrin model, it also predicts orientational diffusion for ellipsoids that |
| 1035 |
< |
exceed explicitly solvated correlation times by a factor of two. |
| 1025 |
> |
The translational diffusion constants from the microcanonical |
| 1026 |
> |
simulations agree well with the predictions of the Perrin model, |
| 1027 |
> |
although the {\it rotational} correlation times are a factor of 2 |
| 1028 |
> |
shorter than expected from hydrodynamic theory. One explanation for |
| 1029 |
> |
the slower rotation of explicitly-solvated ellipsoids is the |
| 1030 |
> |
possibility that solute-solvent collisions happen at both ends of the |
| 1031 |
> |
solute whenever the principal axis of the ellipsoid is turning. In the |
| 1032 |
> |
upper portion of figure \ref{fig:explanation} we sketch a physical |
| 1033 |
> |
picture of this explanation. Since our Langevin integrator is |
| 1034 |
> |
providing nearly quantitative agreement with the Perrin model, it also |
| 1035 |
> |
predicts orientational diffusion for ellipsoids that exceed explicitly |
| 1036 |
> |
solvated correlation times by a factor of two. |
| 1037 |
|
|
| 1038 |
|
\subsection{Rigid dumbbells} |
| 1039 |
|
Perhaps the only {\it composite} rigid body for which analytic |
| 1159 |
|
|
| 1160 |
|
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
| 1161 |
|
used recently as models for lipid |
| 1162 |
< |
molecules.\cite{SunGezelter08,Ayton01} A reference system composed of |
| 1162 |
> |
molecules.\cite{SunX._jp0762020,Ayton01} A reference system composed of |
| 1163 |
|
a single lipid rigid body embedded in a sea of 1929 solvent particles |
| 1164 |
|
was created and run under a microcanonical ensemble. The resulting |
| 1165 |
|
viscosity of this mixture was 0.349 centipoise (as estimated using |
| 1177 |
|
hydrodynamic tensor) are essentially quantitative when compared with |
| 1178 |
|
the explicit solvent simulations for this model system. |
| 1179 |
|
|
| 1180 |
< |
\subsection{Summary} |
| 1181 |
< |
According to our simulations, the Langevin rigid-body integrator we |
| 1182 |
< |
have developed is a reliable way to replace explicit solvent |
| 1183 |
< |
simulations in cases where the detailed solute-solvent interactions do |
| 1184 |
< |
not greatly impact the forces on the solute. In cases where the |
| 1185 |
< |
dielectric screening of the solvent, or specific solute-solvent |
| 1186 |
< |
interactions become important for structural or dynamic features of |
| 1187 |
< |
the solute molecule, this integrator may be less useful. However, for |
| 1188 |
< |
the kinds of coarse-grained modeling that have become popular in |
| 1189 |
< |
recent years, this integrator may prove itself to be quite valuable. |
| 1180 |
> |
\subsection{Summary of comparisons with explicit solvent simulations} |
| 1181 |
> |
The Langevin rigid-body integrator we have developed is a reliable way |
| 1182 |
> |
to replace explicit solvent simulations in cases where the detailed |
| 1183 |
> |
solute-solvent interactions do not greatly impact the behavior of the |
| 1184 |
> |
solute. As such, it has the potential to greatly increase the length |
| 1185 |
> |
and time scales of coarse grained simulations of large solvated |
| 1186 |
> |
molecules. In cases where the dielectric screening of the solvent, or |
| 1187 |
> |
specific solute-solvent interactions become important for structural |
| 1188 |
> |
or dynamic features of the solute molecule, this integrator may be |
| 1189 |
> |
less useful. However, for the kinds of coarse-grained modeling that |
| 1190 |
> |
have become popular in recent years (ellipsoidal side chains, rigid |
| 1191 |
> |
bodies, and molecular-scale models), this integrator may prove itself |
| 1192 |
> |
to be quite valuable. |
| 1193 |
|
|
| 1194 |
|
\begin{figure} |
| 1195 |
|
\centering |
| 1215 |
|
\caption{Translational diffusion constants (D) for the model systems |
| 1216 |
|
calculated using microcanonical simulations (with explicit solvent), |
| 1217 |
|
theoretical predictions, and Langevin simulations (with implicit solvent). |
| 1218 |
< |
Analytical solutions for the exactly-solved hydrodynamics models are |
| 1219 |
< |
from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936} |
| 1218 |
> |
Analytical solutions for the exactly-solved hydrodynamics models are obtained |
| 1219 |
> |
from: Stokes' law (sphere), and Refs. \citen{Perrin1934} and \citen{Perrin1936} |
| 1220 |
|
(ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq} |
| 1221 |
|
(dumbbell). The other model systems have no known analytic solution. |
| 1222 |
< |
All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= |
| 1222 |
> |
All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= |
| 1223 |
|
$10^{-4}$ \AA$^2$ / fs). } |
| 1224 |
|
\begin{tabular}{lccccccc} |
| 1225 |
|
\hline |
| 1272 |
|
|
| 1273 |
|
\section{Application: A rigid-body lipid bilayer} |
| 1274 |
|
|
| 1275 |
< |
The Langevin dynamics integrator was applied to study the formation of |
| 1276 |
< |
corrugated structures emerging from simulations of the coarse grained |
| 1277 |
< |
lipid molecular models presented above. The initial configuration is |
| 1278 |
< |
taken from our molecular dynamics studies on lipid bilayers with |
| 1279 |
< |
lennard-Jones sphere solvents. The solvent molecules were excluded |
| 1280 |
< |
from the system, and the experimental value for the viscosity of water |
| 1281 |
< |
at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects |
| 1282 |
< |
of the solvent. The absence of explicit solvent molecules and the |
| 1283 |
< |
stability of the integrator allowed us to take timesteps of 50 fs. A |
| 1284 |
< |
total simulation run time of 100 ns was sampled. |
| 1285 |
< |
Fig. \ref{fig:bilayer} shows the configuration of the system after 100 |
| 1286 |
< |
ns, and the ripple structure remains stable during the entire |
| 1287 |
< |
trajectory. Compared with using explicit bead-model solvent |
| 1288 |
< |
molecules, the efficiency of the simulation has increased by an order |
| 1289 |
< |
of magnitude. |
| 1275 |
> |
To test the accuracy and efficiency of the new integrator, we applied |
| 1276 |
> |
it to study the formation of corrugated structures emerging from |
| 1277 |
> |
simulations of the coarse grained lipid molecular models presented |
| 1278 |
> |
above. The initial configuration is taken from earlier molecular |
| 1279 |
> |
dynamics studies on lipid bilayers which had used spherical |
| 1280 |
> |
(Lennard-Jones) solvent particles and moderate (480 solvated lipid |
| 1281 |
> |
molecules) system sizes.\cite{SunX._jp0762020} the solvent molecules |
| 1282 |
> |
were excluded from the system and the box was replicated three times |
| 1283 |
> |
in the x- and y- axes (giving a single simulation cell containing 4320 |
| 1284 |
> |
lipids). The experimental value for the viscosity of water at 20C |
| 1285 |
> |
($\eta = 1.00$ cp) was used with the Langevin integrator to mimic the |
| 1286 |
> |
hydrodynamic effects of the solvent. The absence of explicit solvent |
| 1287 |
> |
molecules and the stability of the integrator allowed us to take |
| 1288 |
> |
timesteps of 50 fs. A simulation run time of 30 ns was sampled to |
| 1289 |
> |
calculate structural properties. Fig. \ref{fig:bilayer} shows the |
| 1290 |
> |
configuration of the system after 30 ns. Structural properties of the |
| 1291 |
> |
bilayer (e.g. the head and body $P_2$ order parameters) are nearly |
| 1292 |
> |
identical to those obtained via solvated molecular dynamics. The |
| 1293 |
> |
ripple structure remained stable during the entire trajectory. |
| 1294 |
> |
Compared with using explicit bead-model solvent molecules, the 30 ns |
| 1295 |
> |
trajectory for 4320 lipids with the Langevin integrator is now {\it |
| 1296 |
> |
faster} on the same hardware than the same length trajectory was for |
| 1297 |
> |
the 480-lipid system previously studied. |
| 1298 |
|
|
| 1299 |
|
\begin{figure} |
| 1300 |
|
\centering |
| 1301 |
|
\includegraphics[width=\linewidth]{bilayer} |
| 1302 |
|
\caption[Snapshot of a bilayer of rigid-body models for lipids]{A |
| 1303 |
< |
snapshot of a bilayer composed of rigid-body models for lipid |
| 1303 |
> |
snapshot of a bilayer composed of 4320 rigid-body models for lipid |
| 1304 |
|
molecules evolving using the Langevin integrator described in this |
| 1305 |
|
work.} \label{fig:bilayer} |
| 1306 |
|
\end{figure} |
| 1307 |
|
|
| 1308 |
|
\section{Conclusions} |
| 1309 |
|
|
| 1310 |
< |
We have presented a new Langevin algorithm by incorporating the |
| 1311 |
< |
hydrodynamics properties of arbitrary shaped molecules into an |
| 1312 |
< |
advanced symplectic integration scheme. Further studies in systems |
| 1313 |
< |
involving banana shaped molecules illustrated that the dynamic |
| 1314 |
< |
properties could be preserved by using this new algorithm as an |
| 1315 |
< |
implicit solvent model. |
| 1310 |
> |
We have presented a new algorithm for carrying out Langevin dynamics |
| 1311 |
> |
simulations on complex rigid bodies by incorporating the hydrodynamic |
| 1312 |
> |
resistance tensors for arbitrary shapes into an advanced symplectic |
| 1313 |
> |
integration scheme. The integrator gives quantitative agreement with |
| 1314 |
> |
both analytic and approximate hydrodynamic theories, and works |
| 1315 |
> |
reasonably well at reproducing the solute dynamical properties |
| 1316 |
> |
(diffusion constants, and orientational relaxation times) from |
| 1317 |
> |
explicitly-solvated simulations. For the cases where there are |
| 1318 |
> |
discrepancies between our Langevin integrator and the explicit solvent |
| 1319 |
> |
simulations, two features of molecular simulations help explain the |
| 1320 |
> |
differences. |
| 1321 |
|
|
| 1322 |
+ |
First, the use of ``stick'' boundary conditions for molecular-sized |
| 1323 |
+ |
solutes in a sea of similarly-sized solvent particles may be |
| 1324 |
+ |
problematic. We are certainly not the first group to notice this |
| 1325 |
+ |
difference between hydrodynamic theories and explicitly-solvated |
| 1326 |
+ |
molecular |
| 1327 |
+ |
simulations.\cite{Schmidt:2004fj,Schmidt:2003kx,Ravichandran:1999fk,TANG:1993lr} |
| 1328 |
+ |
The problem becomes particularly noticable in both the translational |
| 1329 |
+ |
diffusion of the spherical particles and the rotational diffusion of |
| 1330 |
+ |
the ellipsoids. In both of these cases it is clear that the |
| 1331 |
+ |
approximations that go into hydrodynamics are the source of the error, |
| 1332 |
+ |
and not the integrator itself. |
| 1333 |
|
|
| 1334 |
+ |
Second, in the case of structures which have substantial surface area |
| 1335 |
+ |
that is inaccessible to solvent particles, the hydrodynamic theories |
| 1336 |
+ |
(and the Langevin integrator) may overestimate the effects of solvent |
| 1337 |
+ |
friction because they overestimate the exposed surface area of the |
| 1338 |
+ |
rigid body. This is particularly noticable in the rotational |
| 1339 |
+ |
diffusion of the dumbbell model. We believe that given a solvent of |
| 1340 |
+ |
known radius, it may be possible to modify the rough shell approach to |
| 1341 |
+ |
place beads on solvent-accessible surface, instead of on the geometric |
| 1342 |
+ |
surface defined by the van der Waals radii of the components of the |
| 1343 |
+ |
rigid body. Further work to confirm the behavior of this new |
| 1344 |
+ |
approximation is ongoing. |
| 1345 |
+ |
|
| 1346 |
|
\section{Acknowledgments} |
| 1347 |
|
Support for this project was provided by the National Science |
| 1348 |
|
Foundation under grant CHE-0134881. T.L. also acknowledges the |
| 1349 |
< |
financial support from center of applied mathematics at University |
| 1350 |
< |
of Notre Dame. |
| 1349 |
> |
financial support from Center of Applied Mathematics at University of |
| 1350 |
> |
Notre Dame. |
| 1351 |
> |
|
| 1352 |
> |
\end{doublespace} |
| 1353 |
|
\newpage |
| 1354 |
|
|
| 1355 |
< |
\bibliographystyle{jcp} |
| 1355 |
> |
\bibliographystyle{jcp2} |
| 1356 |
|
\bibliography{langevin} |
| 1315 |
– |
|
| 1357 |
|
\end{document} |
