| 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} |
| 5 |
< |
\usepackage{amsmath,bm} |
| 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} |
| 20 |
– |
\renewcommand\citemid{\ } % no comma in optional reference note |
| 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{Teng Lin, Xiuquan Sun 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 |
< |
As alternative to Newtonian dynamics, Langevin dynamics, which |
| 66 |
< |
mimics a simple heat bath with stochastic and dissipative forces, |
| 67 |
< |
has been applied in a variety of studies. The stochastic treatment |
| 68 |
< |
of the solvent enables us to carry out substantially longer time |
| 69 |
< |
simulations. Implicit solvent Langevin dynamics simulations of |
| 70 |
< |
met-enkephalin not only outperform explicit solvent simulations for |
| 71 |
< |
computational efficiency, but also agrees very well with explicit |
| 72 |
< |
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 |
| 75 |
< |
instance, the oscillation power spectrum of nanoparticles from |
| 76 |
< |
Langevin dynamics simulation has the same peak frequencies for |
| 77 |
< |
different wave vectors, which recovers the property of magnetic |
| 78 |
< |
excitations in small finite structures.\cite{Berkov2005a} |
| 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 |
| 70 |
> |
dynamic and structural properties obtained from Langevin simulations |
| 71 |
> |
agree quite well with similar properties obtained from explicit |
| 72 |
> |
solvent simulations. |
| 73 |
|
|
| 74 |
< |
%review rigid body dynamics |
| 75 |
< |
Rigid bodies are frequently involved in the modeling of different |
| 76 |
< |
areas, from engineering, physics, to chemistry. For example, |
| 77 |
< |
missiles and vehicle are usually modeled by rigid bodies. The |
| 78 |
< |
movement of the objects in 3D gaming engine or other physics |
| 79 |
< |
simulator is governed by the rigid body dynamics. In molecular |
| 80 |
< |
simulation, rigid body is used to simplify the model in |
| 87 |
< |
protein-protein docking study{\cite{Gray2003}}. |
| 74 |
> |
Recent examples of the usefulness of Langevin simulations include a |
| 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, Liwo and his coworkers suggest that protein folding |
| 79 |
> |
pathways can be explored within a reasonable amount of |
| 80 |
> |
time.\cite{Liwo2005} |
| 81 |
|
|
| 82 |
< |
It is very important to develop stable and efficient methods to |
| 83 |
< |
integrate the equations of motion for orientational degrees of |
| 84 |
< |
freedom. Euler angles are the natural choice to describe the |
| 85 |
< |
rotational degrees of freedom. However, due to $\frac {1}{sin |
| 86 |
< |
\theta}$ singularities, the numerical integration of corresponding |
| 87 |
< |
equations of these motion is very inefficient and inaccurate. |
| 88 |
< |
Although an alternative integrator using multiple sets of Euler |
| 89 |
< |
angles can overcome this difficulty\cite{Barojas1973}, the |
| 90 |
< |
computational penalty and the loss of angular momentum conservation |
| 91 |
< |
still remain. A singularity-free representation utilizing |
| 99 |
< |
quaternions was developed by Evans in 1977.\cite{Evans1977} |
| 100 |
< |
Unfortunately, this approach used a nonseparable Hamiltonian |
| 101 |
< |
resulting from the quaternion representation, which prevented the |
| 102 |
< |
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} |
| 82 |
> |
The stochastic nature of Langevin dynamics also enhances the sampling |
| 83 |
> |
of the system and increases the probability of crossing energy |
| 84 |
> |
barriers.\cite{Cui2003,Banerjee2004} Combining Langevin dynamics with |
| 85 |
> |
Kramers' theory, Klimov and Thirumalai identified free-energy |
| 86 |
> |
barriers by studying the viscosity dependence of the protein folding |
| 87 |
> |
rates.\cite{Klimov1997} In order to account for solvent induced |
| 88 |
> |
interactions missing from the implicit solvent model, Kaya |
| 89 |
> |
incorporated a desolvation free energy barrier into protein |
| 90 |
> |
folding/unfolding studies and discovered a higher free energy barrier |
| 91 |
> |
between the native and denatured states.\cite{HuseyinKaya07012005} |
| 92 |
|
|
| 93 |
< |
A break-through in geometric literature suggests that, in order to |
| 94 |
< |
develop a long-term integration scheme, one should preserve the |
| 95 |
< |
symplectic structure of the propagator. By introducing a conjugate |
| 96 |
< |
momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's |
| 97 |
< |
equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was |
| 98 |
< |
proposed to evolve the Hamiltonian system in a constraint manifold |
| 99 |
< |
by iteratively satisfying the orthogonality constraint $Q^T Q = 1$. |
| 100 |
< |
An alternative method using the quaternion representation was |
| 101 |
< |
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. |
| 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) + 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 \rho$. Here $\eta$ is the viscosity of the |
| 101 |
> |
implicit solvent, and $\rho$ is the hydrodynamic radius of the atom. |
| 102 |
|
|
| 103 |
< |
%review langevin/browninan dynamics for arbitrarily shaped rigid body |
| 104 |
< |
Combining Langevin or Brownian dynamics with rigid body dynamics, |
| 105 |
< |
one can study slow processes in biomolecular systems. Modeling DNA |
| 106 |
< |
as a chain of rigid beads, which are subject to harmonic potentials |
| 107 |
< |
as well as excluded volume potentials, Mielke and his coworkers |
| 108 |
< |
discovered rapid superhelical stress generations from the stochastic |
| 109 |
< |
simulation of twin supercoiling DNA with response to induced |
| 110 |
< |
torques.\cite{Mielke2004} Membrane fusion is another key biological |
| 111 |
< |
process which controls a variety of physiological functions, such as |
| 112 |
< |
release of neurotransmitters \textit{etc}. A typical fusion event |
| 113 |
< |
happens on the time scale of a millisecond, which is impractical to |
| 114 |
< |
study using atomistic models with newtonian mechanics. With the help |
| 115 |
< |
of coarse-grained rigid body model and stochastic dynamics, the |
| 116 |
< |
fusion pathways were explored by many |
| 117 |
< |
researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the |
| 118 |
< |
difficulty of numerical integration of anisotropic rotation, most of |
| 119 |
< |
the rigid body models are simply modeled using spheres, cylinders, |
| 120 |
< |
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 |
| 103 |
> |
The use of rigid substructures,\cite{Chun:2000fj} |
| 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 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 simple rotation evolution scheme consisting of three |
| 123 |
< |
consecutive rotations.\cite{Fernandes2002} Unfortunately, unexpected |
| 124 |
< |
errors and biases are introduced into the system due to the |
| 125 |
< |
arbitrary order of applying the noncommuting rotation |
| 126 |
< |
operators.\cite{Beard2003} Based on the observation the momentum |
| 127 |
< |
relaxation time is much less than the time step, one may ignore the |
| 128 |
< |
inertia in Brownian dynamics. However, the assumption of zero |
| 129 |
< |
average acceleration is not always true for cooperative motion which |
| 130 |
< |
is common in protein motion. An inertial Brownian dynamics (IBD) was |
| 131 |
< |
proposed to address this issue by adding an inertial correction |
| 132 |
< |
term.\cite{Beard2000} As a complement to IBD which has a lower bound |
| 133 |
< |
in time step because of the inertial relaxation time, long-time-step |
| 134 |
< |
inertial dynamics (LTID) can be used to investigate the inertial |
| 135 |
< |
behavior of the polymer segments in low friction |
| 136 |
< |
regime.\cite{Beard2000} LTID can also deal with the rotational |
| 137 |
< |
dynamics for nonskew bodies without translation-rotation coupling by |
| 138 |
< |
separating the translation and rotation motion and taking advantage |
| 139 |
< |
of the analytical solution of hydrodynamics properties. However, |
| 140 |
< |
typical nonskew bodies like cylinders and ellipsoids are inadequate |
| 141 |
< |
to represent most complex macromolecule assemblies. These intricate |
| 142 |
< |
molecules have been represented by a set of beads and their |
| 143 |
< |
hydrodynamic properties can be calculated using variants on the |
| 144 |
< |
standard hydrodynamic interaction tensors. |
| 122 |
> |
introducing a rotational evolution scheme consisting of three |
| 123 |
> |
consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are |
| 124 |
> |
introduced into the system due to the arbitrary order of applying the |
| 125 |
> |
noncommuting rotation operators.\cite{Beard2003} Based on the |
| 126 |
> |
observation the momentum relaxation time is much less than the time |
| 127 |
> |
step, one may ignore the inertia in Brownian dynamics. However, the |
| 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, |
| 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 hydrodynamic |
| 139 |
> |
properties. However, typical nonskew bodies like cylinders and |
| 140 |
> |
ellipsoids are inadequate to represent most complex macromolecular |
| 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 |
< |
The goal of the present work is to develop a Langevin dynamics |
| 147 |
< |
algorithm for arbitrary-shaped rigid particles by integrating the |
| 148 |
< |
accurate estimation of friction tensor from hydrodynamics theory |
| 149 |
< |
into the sophisticated rigid body dynamics algorithms. |
| 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{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 in molecular dynamics simulations using constraints |
| 155 |
> |
rather than rigid body equations of motion. |
| 156 |
|
|
| 157 |
< |
\section{Computational Methods{\label{methodSec}}} |
| 157 |
> |
Euler angles are a natural choice to describe the rotational degrees |
| 158 |
> |
of freedom. However, due to $\frac{1}{\sin \theta}$ singularities, the |
| 159 |
> |
numerical integration of corresponding equations of these motion can |
| 160 |
> |
become inaccurate (and inefficient). Although the use of multiple |
| 161 |
> |
sets of Euler angles can overcome this problem,\cite{Barojas1973} the |
| 162 |
> |
computational penalty and the loss of angular momentum conservation |
| 163 |
> |
remain. A singularity-free representation utilizing quaternions was |
| 164 |
> |
developed by Evans in 1977.\cite{Evans1977} The Evans quaternion |
| 165 |
> |
approach uses a nonseparable Hamiltonian, and this has prevented |
| 166 |
> |
symplectic algorithms from being utilized until very |
| 167 |
> |
recently.\cite{Miller2002} |
| 168 |
|
|
| 169 |
< |
\subsection{\label{introSection:frictionTensor}Friction Tensor} |
| 170 |
< |
Theoretically, the friction kernel can be determined using the |
| 171 |
< |
velocity autocorrelation function. However, this approach becomes |
| 172 |
< |
impractical when the system becomes more and more complicated. |
| 173 |
< |
Instead, various approaches based on hydrodynamics have been |
| 174 |
< |
developed to calculate the friction coefficients. In general, the |
| 175 |
< |
friction tensor $\Xi$ is a $6\times 6$ matrix given by |
| 176 |
< |
\[ |
| 177 |
< |
\Xi = \left( {\begin{array}{*{20}c} |
| 178 |
< |
{\Xi _{}^{tt} } & {\Xi _{}^{rt} } \\ |
| 179 |
< |
{\Xi _{}^{tr} } & {\Xi _{}^{rr} } \\ |
| 180 |
< |
\end{array}} \right). |
| 181 |
< |
\] |
| 182 |
< |
Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$ |
| 183 |
< |
translational friction tensor and rotational resistance (friction) |
| 184 |
< |
tensor respectively, while ${\Xi^{tr} }$ is translation-rotation |
| 185 |
< |
coupling tensor and $ {\Xi^{rt} }$ is rotation-translation coupling |
| 186 |
< |
tensor. When a particle moves in a fluid, it may experience friction |
| 187 |
< |
force or torque along the opposite direction of the velocity or |
| 188 |
< |
angular velocity, |
| 189 |
< |
\[ |
| 169 |
> |
Another approach is the application of holonomic constraints to the |
| 170 |
> |
atoms belonging to the rigid body. Each atom moves independently |
| 171 |
> |
under the normal forces deriving from potential energy and constraints |
| 172 |
> |
are used to guarantee rigidity. However, due to their iterative |
| 173 |
> |
nature, the SHAKE and RATTLE algorithms converge very slowly when the |
| 174 |
> |
number of constraints (and the number of particles that belong to the |
| 175 |
> |
rigid body) increases.\cite{Ryckaert1977,Andersen1983} |
| 176 |
> |
|
| 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 |
> |
${\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 ${\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 |
| 193 |
> |
accurate estimate of the friction tensor from hydrodynamics theory |
| 194 |
> |
into a symplectic rigid body dynamics propagator. In the sections |
| 195 |
> |
below, we review some of the theory of hydrodynamic tensors developed |
| 196 |
> |
primarily for Brownian simulations of multi-particle systems, we then |
| 197 |
> |
present our integration method for a set of generalized Langevin |
| 198 |
> |
equations of motion, and we compare the behavior of the new Langevin |
| 199 |
> |
integrator to dynamical quantities obtained via explicit solvent |
| 200 |
> |
molecular dynamics. |
| 201 |
> |
|
| 202 |
> |
\subsection{\label{introSection:frictionTensor}The Friction Tensor} |
| 203 |
> |
Theoretically, a complete friction kernel for a solute particle can be |
| 204 |
> |
determined using the velocity autocorrelation function from a |
| 205 |
> |
simulation with explicit solvent molecules. However, this approach |
| 206 |
> |
becomes impractical when the solute becomes complex. Instead, various |
| 207 |
> |
approaches based on hydrodynamics have been developed to calculate |
| 208 |
> |
static friction coefficients. In general, the friction tensor $\Xi$ is |
| 209 |
> |
a $6\times 6$ matrix given by |
| 210 |
> |
\begin{equation} |
| 211 |
> |
\Xi = \left( \begin{array}{*{20}c} |
| 212 |
> |
\Xi^{tt} & \Xi^{rt} \\ |
| 213 |
> |
\Xi^{tr} & \Xi^{rr} \\ |
| 214 |
> |
\end{array} \right). |
| 215 |
> |
\end{equation} |
| 216 |
> |
Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and |
| 217 |
> |
rotational resistance (friction) tensors respectively, while |
| 218 |
> |
$\Xi^{tr}$ is translation-rotation coupling tensor and $\Xi^{rt}$ is |
| 219 |
> |
rotation-translation coupling tensor. When a particle moves in a |
| 220 |
> |
fluid, it may experience a friction force ($\mathbf{f}_f$) and torque |
| 221 |
> |
($\mathbf{\tau}_f$) in opposition to the velocity ($\mathbf{v}$) and |
| 222 |
> |
body-fixed angular velocity ($\mathbf{\omega}$), |
| 223 |
> |
\begin{equation} |
| 224 |
|
\left( \begin{array}{l} |
| 225 |
< |
F_R \\ |
| 226 |
< |
\tau _R \\ |
| 227 |
< |
\end{array} \right) = - \left( {\begin{array}{*{20}c} |
| 228 |
< |
{\Xi ^{tt} } & {\Xi ^{rt} } \\ |
| 229 |
< |
{\Xi ^{tr} } & {\Xi ^{rr} } \\ |
| 230 |
< |
\end{array}} \right)\left( \begin{array}{l} |
| 231 |
< |
v \\ |
| 232 |
< |
w \\ |
| 233 |
< |
\end{array} \right) |
| 234 |
< |
\] |
| 208 |
< |
where $F_r$ is the friction force and $\tau _R$ is the friction |
| 209 |
< |
torque. |
| 225 |
> |
\mathbf{f}_f \\ |
| 226 |
> |
\mathbf{\tau}_f \\ |
| 227 |
> |
\end{array} \right) = - \left( \begin{array}{*{20}c} |
| 228 |
> |
\Xi^{tt} & \Xi^{rt} \\ |
| 229 |
> |
\Xi^{tr} & \Xi^{rr} \\ |
| 230 |
> |
\end{array} \right)\left( \begin{array}{l} |
| 231 |
> |
\mathbf{v} \\ |
| 232 |
> |
\mathbf{\omega} \\ |
| 233 |
> |
\end{array} \right). |
| 234 |
> |
\end{equation} |
| 235 |
|
|
| 236 |
|
\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} |
| 237 |
< |
|
| 238 |
< |
For a spherical particle with slip boundary conditions, the |
| 239 |
< |
translational and rotational friction constant can be calculated |
| 240 |
< |
from Stoke's law, |
| 241 |
< |
\[ |
| 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} \\ |
| 246 |
< |
\end{array}} \right) |
| 247 |
< |
\] |
| 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 \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 |
< |
\[ |
| 250 |
< |
\Xi ^{rr} = \left( {\begin{array}{*{20}c} |
| 251 |
< |
{8\pi \eta R^3 } & 0 & 0 \\ |
| 252 |
< |
0 & {8\pi \eta R^3 } & 0 \\ |
| 253 |
< |
0 & 0 & {8\pi \eta R^3 } \\ |
| 254 |
< |
\end{array}} \right) |
| 255 |
< |
\] |
| 256 |
< |
where $\eta$ is the viscosity of the solvent and $R$ is the |
| 257 |
< |
hydrodynamic radius. |
| 249 |
> |
\begin{equation} |
| 250 |
> |
\label{eq:StokesRotation} |
| 251 |
> |
\Xi^{rr} = \left( \begin{array}{*{20}c} |
| 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 $\rho$ is the |
| 258 |
> |
hydrodynamic radius. The presence of the rotational resistance tensor |
| 259 |
> |
implies that the spherical body has internal structure and |
| 260 |
> |
orientational degrees of freedom that must be propagated in time. For |
| 261 |
> |
non-structured spherical bodies (i.e. the atoms in a traditional |
| 262 |
> |
molecular dynamics simulation) these degrees of freedom do not exist. |
| 263 |
|
|
| 264 |
|
Other non-spherical shapes, such as cylinders and ellipsoids, are |
| 265 |
< |
widely used as references for developing new hydrodynamics theory, |
| 265 |
> |
widely used as references for developing new hydrodynamic theories, |
| 266 |
|
because their properties can be calculated exactly. In 1936, Perrin |
| 267 |
< |
extended Stokes's law to general ellipsoids, also called a triaxial |
| 268 |
< |
ellipsoid, which is given in Cartesian coordinates |
| 269 |
< |
by\cite{Perrin1934, Perrin1936} |
| 270 |
< |
\[ |
| 271 |
< |
\frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2 |
| 272 |
< |
}} = 1 |
| 273 |
< |
\] |
| 274 |
< |
where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately, |
| 275 |
< |
due to the complexity of the elliptic integral, only the ellipsoid |
| 276 |
< |
with the restriction of two axes being equal, \textit{i.e.} |
| 277 |
< |
prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved |
| 278 |
< |
exactly. Introducing an elliptic integral parameter $S$ for prolate |
| 279 |
< |
ellipsoids : |
| 280 |
< |
\[ |
| 281 |
< |
S = \frac{2}{{\sqrt {a^2 - b^2 } }}\ln \frac{{a + \sqrt {a^2 - b^2 |
| 282 |
< |
} }}{b}, |
| 283 |
< |
\] |
| 284 |
< |
and oblate ellipsoids: |
| 285 |
< |
\[ |
| 286 |
< |
S = \frac{2}{{\sqrt {b^2 - a^2 } }}arctg\frac{{\sqrt {b^2 - a^2 } |
| 287 |
< |
}}{a}, |
| 288 |
< |
\] |
| 289 |
< |
one can write down the translational and rotational resistance |
| 290 |
< |
tensors |
| 291 |
< |
\begin{eqnarray*} |
| 292 |
< |
\Xi _a^{tt} & = & 16\pi \eta \frac{{a^2 - b^2 }}{{(2a^2 - b^2 )S - 2a}}. \\ |
| 293 |
< |
\Xi _b^{tt} & = & \Xi _c^{tt} = 32\pi \eta \frac{{a^2 - b^2 }}{{(2a^2 - 3b^2 )S + |
| 294 |
< |
2a}}, |
| 295 |
< |
\end{eqnarray*} |
| 296 |
< |
and |
| 297 |
< |
\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}}. |
| 270 |
< |
\end{eqnarray*} |
| 267 |
> |
extended Stokes' law to general |
| 268 |
> |
ellipsoids,\cite{Perrin1934,Perrin1936} described in Cartesian |
| 269 |
> |
coordinates as |
| 270 |
> |
\begin{equation} |
| 271 |
> |
\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1. |
| 272 |
> |
\end{equation} |
| 273 |
> |
Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the |
| 274 |
> |
complexity of the elliptic integral, only uniaxial ellipsoids, either |
| 275 |
> |
prolate ($a \ge b = c$) or oblate ($a < b = c$), were solved |
| 276 |
> |
exactly. Introducing an elliptic integral parameter $S$ for prolate, |
| 277 |
> |
\begin{equation} |
| 278 |
> |
S = \frac{2}{\sqrt{a^2 - b^2}} \ln \frac{a + \sqrt{a^2 - b^2}}{b}, |
| 279 |
> |
\end{equation} |
| 280 |
> |
and oblate, |
| 281 |
> |
\begin{equation} |
| 282 |
> |
S = \frac{2}{\sqrt {b^2 - a^2 }} \arctan \frac{\sqrt {b^2 - a^2}}{a}, |
| 283 |
> |
\end{equation} |
| 284 |
> |
ellipsoids, it is possible to write down exact solutions for the |
| 285 |
> |
resistance tensors. As is the case for spherical bodies, the translational, |
| 286 |
> |
\begin{eqnarray} |
| 287 |
> |
\Xi_a^{tt} & = & 16\pi \eta \frac{a^2 - b^2}{(2a^2 - b^2 )S - 2a}. \\ |
| 288 |
> |
\Xi_b^{tt} = \Xi_c^{tt} & = & 32\pi \eta \frac{a^2 - b^2 }{(2a^2 - 3b^2 )S + 2a}, |
| 289 |
> |
\end{eqnarray} |
| 290 |
> |
and rotational, |
| 291 |
> |
\begin{eqnarray} |
| 292 |
> |
\Xi_a^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^2 - b^2 )b^2}{2a - b^2 S}, \\ |
| 293 |
> |
\Xi_b^{rr} = \Xi_c^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^4 - b^4)}{(2a^2 - b^2 )S - 2a} |
| 294 |
> |
\end{eqnarray} |
| 295 |
> |
resistance tensors are diagonal $3 \times 3$ matrices. For both |
| 296 |
> |
spherical and ellipsoidal particles, the translation-rotation and |
| 297 |
> |
rotation-translation coupling tensors are zero. |
| 298 |
|
|
| 299 |
|
\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}} |
| 300 |
+ |
There is no analytical solution for the friction tensor for rigid |
| 301 |
+ |
molecules of arbitrary shape. The ellipsoid of revolution and general |
| 302 |
+ |
triaxial ellipsoid models have been widely used to approximate the |
| 303 |
+ |
hydrodynamic properties of rigid bodies. However, the mapping from all |
| 304 |
+ |
possible ellipsoidal spaces ($r$-space) to all possible combinations |
| 305 |
+ |
of rotational diffusion coefficients ($D$-space) is not |
| 306 |
+ |
unique.\cite{Wegener1979} Additionally, because there is intrinsic |
| 307 |
+ |
coupling between translational and rotational motion of {\it skew} |
| 308 |
+ |
rigid bodies, general ellipsoids are not always suitable for modeling |
| 309 |
+ |
rigid molecules. A number of studies have been devoted to determining |
| 310 |
+ |
the friction tensor for irregular shapes using methods in which the |
| 311 |
+ |
molecule of interest is modeled with a combination of |
| 312 |
+ |
spheres\cite{Carrasco1999} and the hydrodynamic properties of the |
| 313 |
+ |
molecule are then calculated using a set of two-point interaction |
| 314 |
+ |
tensors. We have found the {\it bead} and {\it rough shell} models of |
| 315 |
+ |
Carrasco and Garc\'{i}a de la Torre to be the most useful of these |
| 316 |
+ |
methods,\cite{Carrasco1999} and we review the basic outline of the |
| 317 |
+ |
rough shell approach here. A more thorough explanation can be found |
| 318 |
+ |
in Ref. \citen{Carrasco1999}. |
| 319 |
|
|
| 320 |
< |
Unlike spherical and other simply shaped molecules, there is no |
| 321 |
< |
analytical solution for the friction tensor for arbitrarily shaped |
| 322 |
< |
rigid molecules. The ellipsoid of revolution model and general |
| 323 |
< |
triaxial ellipsoid model have been used to approximate the |
| 324 |
< |
hydrodynamic properties of rigid bodies. However, since the mapping |
| 279 |
< |
from all possible ellipsoidal spaces, $r$-space, to all possible |
| 280 |
< |
combination of rotational diffusion coefficients, $D$-space, is not |
| 281 |
< |
unique\cite{Wegener1979} as well as the intrinsic coupling between |
| 282 |
< |
translational and rotational motion of rigid bodies, general |
| 283 |
< |
ellipsoids are not always suitable for modeling arbitrarily shaped |
| 284 |
< |
rigid molecules. A number of studies have been devoted to |
| 285 |
< |
determining the friction tensor for irregularly shaped rigid bodies |
| 286 |
< |
using more advanced methods where the molecule of interest was |
| 287 |
< |
modeled by a combinations of spheres\cite{Carrasco1999} and the |
| 288 |
< |
hydrodynamics properties of the molecule can be calculated using the |
| 289 |
< |
hydrodynamic interaction tensor. Let us consider a rigid assembly of |
| 290 |
< |
$N$ beads immersed in a continuous medium. Due to hydrodynamic |
| 291 |
< |
interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different |
| 292 |
< |
than its unperturbed velocity $v_i$, |
| 293 |
< |
\[ |
| 294 |
< |
v'_i = v_i - \sum\limits_{j \ne i} {T_{ij} F_j } |
| 295 |
< |
\] |
| 296 |
< |
where $F_i$ is the frictional force, and $T_{ij}$ is the |
| 297 |
< |
hydrodynamic interaction tensor. The friction force of $i$th bead is |
| 298 |
< |
proportional to its ``net'' velocity |
| 320 |
> |
Consider a rigid assembly of $N$ small beads moving through a |
| 321 |
> |
continuous medium. Due to hydrodynamic interactions between the |
| 322 |
> |
beads, the net velocity of the $i^\mathrm{th}$ bead relative to the |
| 323 |
> |
medium, ${\bf v}'_i$, is different than its unperturbed velocity ${\bf |
| 324 |
> |
v}_i$, |
| 325 |
|
\begin{equation} |
| 326 |
< |
F_i = \zeta _i v_i - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }. |
| 301 |
< |
\label{introEquation:tensorExpression} |
| 326 |
> |
{\bf v}'_i = {\bf v}_i - \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j } |
| 327 |
|
\end{equation} |
| 328 |
< |
This equation is the basis for deriving the hydrodynamic tensor. In |
| 329 |
< |
1930, Oseen and Burgers gave a simple solution to |
| 330 |
< |
Eq.~\ref{introEquation:tensorExpression} |
| 328 |
> |
where ${\bf F}_j$ is the frictional force on the medium due to bead $j$, and |
| 329 |
> |
${\bf T}_{ij}$ is the hydrodynamic interaction tensor between the two beads. |
| 330 |
> |
The frictional force felt by the $i^\mathrm{th}$ bead is proportional to |
| 331 |
> |
its net velocity |
| 332 |
|
\begin{equation} |
| 333 |
< |
T_{ij} = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij} |
| 334 |
< |
R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor} |
| 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 |
< |
Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$. |
| 337 |
< |
A second order expression for element of different size was |
| 338 |
< |
introduced by Rotne and Prager\cite{Rotne1969} and improved by |
| 336 |
> |
Eq. (\ref{introEquation:tensorExpression}) defines the two-point |
| 337 |
> |
hydrodynamic tensor, ${\bf T}_{ij}$. There have been many proposed |
| 338 |
> |
solutions to this equation, including the simple solution given by |
| 339 |
> |
Oseen and Burgers in 1930 for two beads of identical radius. A second |
| 340 |
> |
order expression for beads of different hydrodynamic radii was |
| 341 |
> |
introduced by Rotne and Prager,\cite{Rotne1969} and improved by |
| 342 |
|
Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977} |
| 343 |
|
\begin{equation} |
| 344 |
< |
T_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I + |
| 345 |
< |
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma |
| 346 |
< |
_i^2 + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} - |
| 347 |
< |
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. |
| 344 |
> |
{\bf T}_{ij} = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {{\bf I} + |
| 345 |
> |
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\rho |
| 346 |
> |
_i^2 + \rho_j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} - |
| 347 |
> |
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. |
| 348 |
|
\label{introEquation:RPTensorNonOverlapped} |
| 349 |
|
\end{equation} |
| 350 |
< |
Both of the Eq.~\ref{introEquation:oseenTensor} and |
| 351 |
< |
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption |
| 352 |
< |
$R_{ij} \ge \sigma _i + \sigma _j$. An alternative expression for |
| 353 |
< |
overlapping beads with the same radius, $\sigma$, is given by |
| 350 |
> |
Here ${\bf R}_{ij}$ is the distance vector between beads $i$ and $j$. Both |
| 351 |
> |
the Oseen-Burgers tensor and |
| 352 |
> |
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption that |
| 353 |
> |
the beads do not overlap ($R_{ij} \ge \rho_i + \rho_j$). |
| 354 |
> |
|
| 355 |
> |
To calculate the resistance tensor for a body represented as the union |
| 356 |
> |
of many non-overlapping beads, we first pick an arbitrary origin $O$ |
| 357 |
> |
and then construct a $3N \times 3N$ supermatrix consisting of $N |
| 358 |
> |
\times N$ ${\bf B}_{ij}$ blocks |
| 359 |
|
\begin{equation} |
| 360 |
< |
T_{ij} = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 - |
| 361 |
< |
\frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I + |
| 362 |
< |
\frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right] |
| 363 |
< |
\label{introEquation:RPTensorOverlapped} |
| 360 |
> |
{\bf B} = \left( \begin{array}{*{20}c} |
| 361 |
> |
{\bf B}_{11} & \ldots & {\bf B}_{1N} \\ |
| 362 |
> |
\vdots & \ddots & \vdots \\ |
| 363 |
> |
{\bf B}_{N1} & \cdots & {\bf B}_{NN} |
| 364 |
> |
\end{array} \right) |
| 365 |
|
\end{equation} |
| 366 |
< |
To calculate the resistance tensor at an arbitrary origin $O$, we |
| 367 |
< |
construct a $3N \times 3N$ matrix consisting of $N \times N$ |
| 333 |
< |
$B_{ij}$ blocks |
| 366 |
> |
${\bf B}_{ij}$ is a version of the hydrodynamic tensor which includes the |
| 367 |
> |
self-contributions for spheres, |
| 368 |
|
\begin{equation} |
| 369 |
< |
B = \left( {\begin{array}{*{20}c} |
| 370 |
< |
{B_{11} } & \ldots & {B_{1N} } \\ |
| 337 |
< |
\vdots & \ddots & \vdots \\ |
| 338 |
< |
{B_{N1} } & \cdots & {B_{NN} } \\ |
| 339 |
< |
\end{array}} \right), |
| 369 |
> |
{\bf B}_{ij} = \delta _{ij} \frac{{\bf I}}{{6\pi \eta R_{ij}}} + (1 - \delta_{ij} |
| 370 |
> |
){\bf T}_{ij} |
| 371 |
|
\end{equation} |
| 372 |
< |
where $B_{ij}$ is given by |
| 373 |
< |
\[ |
| 374 |
< |
B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij} |
| 375 |
< |
)T_{ij} |
| 376 |
< |
\] |
| 377 |
< |
where $\delta _{ij}$ is the Kronecker delta function. Inverting the |
| 378 |
< |
$B$ matrix, we obtain |
| 379 |
< |
\[ |
| 380 |
< |
C = B^{ - 1} = \left( {\begin{array}{*{20}c} |
| 381 |
< |
{C_{11} } & \ldots & {C_{1N} } \\ |
| 382 |
< |
\vdots & \ddots & \vdots \\ |
| 383 |
< |
{C_{N1} } & \cdots & {C_{NN} } \\ |
| 384 |
< |
\end{array}} \right), |
| 385 |
< |
\] |
| 386 |
< |
which can be partitioned into $N \times N$ $3 \times 3$ block |
| 387 |
< |
$C_{ij}$. With the help of $C_{ij}$ and the skew matrix $U_i$ |
| 388 |
< |
\[ |
| 389 |
< |
U_i = \left( {\begin{array}{*{20}c} |
| 390 |
< |
0 & { - z_i } & {y_i } \\ |
| 391 |
< |
{z_i } & 0 & { - x_i } \\ |
| 361 |
< |
{ - y_i } & {x_i } & 0 \\ |
| 362 |
< |
\end{array}} \right) |
| 363 |
< |
\] |
| 372 |
> |
where $\delta_{ij}$ is the Kronecker delta function. Inverting the |
| 373 |
> |
${\bf B}$ matrix, we obtain |
| 374 |
> |
\begin{equation} |
| 375 |
> |
{\bf C} = {\bf B}^{ - 1} = \left(\begin{array}{*{20}c} |
| 376 |
> |
{\bf C}_{11} & \ldots & {\bf C}_{1N} \\ |
| 377 |
> |
\vdots & \ddots & \vdots \\ |
| 378 |
> |
{\bf C}_{N1} & \cdots & {\bf C}_{NN} |
| 379 |
> |
\end{array} \right), |
| 380 |
> |
\end{equation} |
| 381 |
> |
which can be partitioned into $N \times N$ blocks labeled ${\bf C}_{ij}$. |
| 382 |
> |
(Each of the ${\bf C}_{ij}$ blocks is a $3 \times 3$ matrix.) Using the |
| 383 |
> |
skew matrix, |
| 384 |
> |
\begin{equation} |
| 385 |
> |
{\bf U}_i = \left(\begin{array}{*{20}c} |
| 386 |
> |
0 & -z_i & y_i \\ |
| 387 |
> |
z_i & 0 & - x_i \\ |
| 388 |
> |
-y_i & x_i & 0 |
| 389 |
> |
\end{array}\right) |
| 390 |
> |
\label{eq:skewMatrix} |
| 391 |
> |
\end{equation} |
| 392 |
|
where $x_i$, $y_i$, $z_i$ are the components of the vector joining |
| 393 |
< |
bead $i$ and origin $O$, the elements of resistance tensor at |
| 394 |
< |
arbitrary origin $O$ can be written as |
| 393 |
> |
bead $i$ and origin $O$, the elements of the resistance tensor (at the |
| 394 |
> |
arbitrary origin $O$) can be written as |
| 395 |
|
\begin{eqnarray} |
| 368 |
– |
\Xi _{}^{tt} & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\ |
| 369 |
– |
\Xi _{}^{tr} & = & \Xi _{}^{rt} = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\ |
| 370 |
– |
\Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\ |
| 396 |
|
\label{introEquation:ResistanceTensorArbitraryOrigin} |
| 397 |
+ |
\Xi^{tt} & = & \sum\limits_i {\sum\limits_j {{\bf C}_{ij} } } \notag , \\ |
| 398 |
+ |
\Xi^{tr} = \Xi _{}^{rt} & = & \sum\limits_i {\sum\limits_j {{\bf U}_i {\bf C}_{ij} } } , \\ |
| 399 |
+ |
\Xi^{rr} & = & -\sum\limits_i \sum\limits_j {\bf U}_i {\bf C}_{ij} {\bf U}_j + 6 \eta V {\bf I}. \notag |
| 400 |
|
\end{eqnarray} |
| 401 |
< |
The resistance tensor depends on the origin to which they refer. The |
| 402 |
< |
proper location for applying the friction force is the center of |
| 403 |
< |
resistance (or center of reaction), at which the trace of rotational |
| 404 |
< |
resistance tensor, $ \Xi ^{rr}$ reaches a minimum value. |
| 377 |
< |
Mathematically, the center of resistance is defined as an unique |
| 378 |
< |
point of the rigid body at which the translation-rotation coupling |
| 379 |
< |
tensors are symmetric, |
| 401 |
> |
The final term in the expression for $\Xi^{rr}$ is a correction that |
| 402 |
> |
accounts for errors in the rotational motion of the bead models. The |
| 403 |
> |
additive correction uses the solvent viscosity ($\eta$) as well as the |
| 404 |
> |
total volume of the beads that contribute to the hydrodynamic model, |
| 405 |
|
\begin{equation} |
| 406 |
< |
\Xi^{tr} = \left( {\Xi^{tr} } \right)^T |
| 382 |
< |
\label{introEquation:definitionCR} |
| 406 |
> |
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3, |
| 407 |
|
\end{equation} |
| 408 |
< |
From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, |
| 409 |
< |
we can easily derive that the translational resistance tensor is |
| 410 |
< |
origin independent, while the rotational resistance tensor and |
| 411 |
< |
translation-rotation coupling resistance tensor depend on the |
| 412 |
< |
origin. Given the resistance tensor at an arbitrary origin $O$, and |
| 413 |
< |
a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can |
| 414 |
< |
obtain the resistance tensor at $P$ by |
| 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} |
| 412 |
> |
|
| 413 |
> |
In general, resistance tensors depend on the origin at which they were |
| 414 |
> |
computed. However, the proper location for applying the friction |
| 415 |
> |
force is the center of resistance, the special point at which the |
| 416 |
> |
trace of rotational resistance tensor, $\Xi^{rr}$ reaches a minimum |
| 417 |
> |
value. Mathematically, the center of resistance can also be defined |
| 418 |
> |
as the unique point for a rigid body at which the translation-rotation |
| 419 |
> |
coupling tensors are symmetric, |
| 420 |
|
\begin{equation} |
| 421 |
< |
\begin{array}{l} |
| 422 |
< |
\Xi _P^{tt} = \Xi _O^{tt} \\ |
| 394 |
< |
\Xi _P^{tr} = \Xi _P^{rt} = \Xi _O^{tr} - U_{OP} \Xi _O^{tt} \\ |
| 395 |
< |
\Xi _P^{rr} = \Xi _O^{rr} - U_{OP} \Xi _O^{tt} U_{OP} + \Xi _O^{tr} U_{OP} - U_{OP} \Xi _O^{{tr} ^{^T }} \\ |
| 396 |
< |
\end{array} |
| 397 |
< |
\label{introEquation:resistanceTensorTransformation} |
| 421 |
> |
\Xi^{tr} = \left(\Xi^{tr}\right)^T |
| 422 |
> |
\label{introEquation:definitionCR} |
| 423 |
|
\end{equation} |
| 424 |
< |
where |
| 425 |
< |
\[ |
| 426 |
< |
U_{OP} = \left( {\begin{array}{*{20}c} |
| 427 |
< |
0 & { - z_{OP} } & {y_{OP} } \\ |
| 428 |
< |
{z_i } & 0 & { - x_{OP} } \\ |
| 429 |
< |
{ - y_{OP} } & {x_{OP} } & 0 \\ |
| 430 |
< |
\end{array}} \right) |
| 431 |
< |
\] |
| 432 |
< |
Using Eq.~\ref{introEquation:definitionCR} and |
| 433 |
< |
Eq.~\ref{introEquation:resistanceTensorTransformation}, one can |
| 434 |
< |
locate the position of center of resistance, |
| 435 |
< |
\begin{eqnarray*} |
| 436 |
< |
\left( \begin{array}{l} |
| 437 |
< |
x_{OR} \\ |
| 438 |
< |
y_{OR} \\ |
| 439 |
< |
z_{OR} \\ |
| 440 |
< |
\end{array} \right) & = &\left( {\begin{array}{*{20}c} |
| 441 |
< |
{(\Xi _O^{rr} )_{yy} + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} } \\ |
| 442 |
< |
{ - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz} + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} } \\ |
| 443 |
< |
{ - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx} + (\Xi _O^{rr} )_{yy} } \\ |
| 444 |
< |
\end{array}} \right)^{ - 1} \\ |
| 445 |
< |
& & \left( \begin{array}{l} |
| 446 |
< |
(\Xi _O^{tr} )_{yz} - (\Xi _O^{tr} )_{zy} \\ |
| 447 |
< |
(\Xi _O^{tr} )_{zx} - (\Xi _O^{tr} )_{xz} \\ |
| 448 |
< |
(\Xi _O^{tr} )_{xy} - (\Xi _O^{tr} )_{yx} \\ |
| 449 |
< |
\end{array} \right) \\ |
| 450 |
< |
\end{eqnarray*} |
| 451 |
< |
where $x_OR$, $y_OR$, $z_OR$ are the components of the vector |
| 424 |
> |
From Eq. \ref{introEquation:ResistanceTensorArbitraryOrigin}, we can |
| 425 |
> |
easily derive that the {\it translational} resistance tensor is origin |
| 426 |
> |
independent, while the rotational resistance tensor and |
| 427 |
> |
translation-rotation coupling resistance tensor depend on the |
| 428 |
> |
origin. Given the resistance tensor at an arbitrary origin $O$, and a |
| 429 |
> |
vector ,${\bf r}_{OP} = (x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we |
| 430 |
> |
can obtain the resistance tensor at $P$ by |
| 431 |
> |
\begin{eqnarray} |
| 432 |
> |
\label{introEquation:resistanceTensorTransformation} |
| 433 |
> |
\Xi_P^{tt} & = & \Xi_O^{tt} \notag \\ |
| 434 |
> |
\Xi_P^{tr} = \Xi_P^{rt} & = & \Xi_O^{tr} - {\bf U}_{OP} \Xi _O^{tt} \\ |
| 435 |
> |
\Xi_P^{rr} & = &\Xi_O^{rr} - {\bf U}_{OP} \Xi_O^{tt} {\bf U}_{OP} |
| 436 |
> |
+ \Xi_O^{tr} {\bf U}_{OP} - {\bf U}_{OP} \left( \Xi_O^{tr} |
| 437 |
> |
\right)^{^T} \notag |
| 438 |
> |
\end{eqnarray} |
| 439 |
> |
where ${\bf U}_{OP}$ is the skew matrix (Eq. (\ref{eq:skewMatrix})) |
| 440 |
> |
for the vector between the origin $O$ and the point $P$. Using |
| 441 |
> |
Eqs.~\ref{introEquation:definitionCR}~and~\ref{introEquation:resistanceTensorTransformation}, |
| 442 |
> |
one can locate the position of center of resistance, |
| 443 |
> |
\begin{equation*} |
| 444 |
> |
\left(\begin{array}{l} |
| 445 |
> |
x_{OR} \\ |
| 446 |
> |
y_{OR} \\ |
| 447 |
> |
z_{OR} |
| 448 |
> |
\end{array}\right) = |
| 449 |
> |
\left(\begin{array}{*{20}c} |
| 450 |
> |
(\Xi_O^{rr})_{yy} + (\Xi_O^{rr})_{zz} & -(\Xi_O^{rr})_{xy} & -(\Xi_O^{rr})_{xz} \\ |
| 451 |
> |
-(\Xi_O^{rr})_{xy} & (\Xi_O^{rr})_{zz} + (\Xi_O^{rr})_{xx} & -(\Xi_O^{rr})_{yz} \\ |
| 452 |
> |
-(\Xi_O^{rr})_{xz} & -(\Xi_O^{rr})_{yz} & (\Xi_O^{rr})_{xx} + (\Xi_O^{rr})_{yy} \\ |
| 453 |
> |
\end{array}\right)^{-1} |
| 454 |
> |
\left(\begin{array}{l} |
| 455 |
> |
(\Xi_O^{tr})_{yz} - (\Xi_O^{tr})_{zy} \\ |
| 456 |
> |
(\Xi_O^{tr})_{zx} - (\Xi_O^{tr})_{xz} \\ |
| 457 |
> |
(\Xi_O^{tr})_{xy} - (\Xi_O^{tr})_{yx} |
| 458 |
> |
\end{array}\right) |
| 459 |
> |
\end{equation*} |
| 460 |
> |
where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector |
| 461 |
|
joining center of resistance $R$ and origin $O$. |
| 462 |
|
|
| 463 |
< |
\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
| 463 |
> |
For a general rigid molecular substructure, finding the $6 \times 6$ |
| 464 |
> |
resistance tensor can be a computationally demanding task. First, a |
| 465 |
> |
lattice of small beads that extends well beyond the boundaries of the |
| 466 |
> |
rigid substructure is created. The lattice is typically composed of |
| 467 |
> |
0.25 \AA\ beads on a dense FCC lattice. The lattice constant is taken |
| 468 |
> |
to be the bead diameter, so that adjacent beads are touching, but do |
| 469 |
> |
not overlap. To make a shape corresponding to the rigid structure, |
| 470 |
> |
beads that sit on lattice sites that are outside the van der Waals |
| 471 |
> |
radii of all of the atoms comprising the rigid body are excluded from |
| 472 |
> |
the calculation. |
| 473 |
|
|
| 474 |
+ |
For large structures, most of the beads will be deep within the rigid |
| 475 |
+ |
body and will not contribute to the hydrodynamic tensor. In the {\it |
| 476 |
+ |
rough shell} approach, beads which have all of their lattice neighbors |
| 477 |
+ |
inside the structure are considered interior beads, and are removed |
| 478 |
+ |
from the calculation. After following this procedure, only those |
| 479 |
+ |
beads in direct contact with the van der Waals surface of the rigid |
| 480 |
+ |
body are retained. For reasonably large molecular structures, this |
| 481 |
+ |
truncation can still produce bead assemblies with thousands of |
| 482 |
+ |
members. |
| 483 |
+ |
|
| 484 |
+ |
If all of the {\it atoms} comprising the rigid substructure are |
| 485 |
+ |
spherical and non-overlapping, the tensor in |
| 486 |
+ |
Eq.~(\ref{introEquation:RPTensorNonOverlapped}) may be used directly |
| 487 |
+ |
using the atoms themselves as the hydrodynamic beads. This is a |
| 488 |
+ |
variant of the {\it bead model} approach of Carrasco and Garc\'{i}a de |
| 489 |
+ |
la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$ |
| 490 |
+ |
matrix can be quite small, and the calculation of the hydrodynamic |
| 491 |
+ |
tensor is straightforward. |
| 492 |
+ |
|
| 493 |
+ |
In general, the inversion of the ${\bf B}$ matrix is the most |
| 494 |
+ |
computationally demanding task. This inversion is done only once for |
| 495 |
+ |
each type of rigid structure. We have used straightforward |
| 496 |
+ |
LU-decomposition to solve the linear system and to obtain the elements |
| 497 |
+ |
of ${\bf C}$. Once ${\bf C}$ has been obtained, the location of the |
| 498 |
+ |
center of resistance ($R$) is found and the resistance tensor at this |
| 499 |
+ |
point is calculated. The $3 \times 1$ vector giving the location of |
| 500 |
+ |
the rigid body's center of resistance and the $6 \times 6$ resistance |
| 501 |
+ |
tensor are then stored for use in the Langevin dynamics calculation. |
| 502 |
+ |
These quantities depend on solvent viscosity and temperature and must |
| 503 |
+ |
be recomputed if different simulation conditions are required. |
| 504 |
+ |
|
| 505 |
+ |
\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
| 506 |
+ |
|
| 507 |
|
Consider the Langevin equations of motion in generalized coordinates |
| 508 |
|
\begin{equation} |
| 509 |
< |
M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)} + F_{r,i} (t) |
| 509 |
> |
{\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) + |
| 510 |
> |
{\bf F}_{f}(t) + {\bf F}_{r}(t) |
| 511 |
|
\label{LDGeneralizedForm} |
| 512 |
|
\end{equation} |
| 513 |
< |
where $M_i$ is a $6\times6$ generalized diagonal mass (include mass |
| 514 |
< |
and moment of inertial) matrix and $V_i$ is a generalized velocity, |
| 515 |
< |
$V_i = V_i(v_i,\omega _i)$. The right side of |
| 516 |
< |
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in |
| 517 |
< |
lab-fixed frame, systematic force $F_{s,i}$, dissipative force |
| 518 |
< |
$F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the |
| 519 |
< |
system in Newtownian mechanics typically refers to lab-fixed frame, |
| 520 |
< |
it is also convenient to handle the rotation of rigid body in |
| 521 |
< |
body-fixed frame. Thus the friction and random forces are calculated |
| 522 |
< |
in body-fixed frame and converted back to lab-fixed frame by: |
| 523 |
< |
\[ |
| 524 |
< |
\begin{array}{l} |
| 525 |
< |
F_{f,i}^l (t) = Q^T F_{f,i}^b (t), \\ |
| 526 |
< |
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$ |
| 513 |
> |
where ${\bf M}$ is a $6 \times 6$ diagonal mass matrix (which |
| 514 |
> |
includes the mass of the rigid body as well as the moments of inertia |
| 515 |
> |
in the body-fixed frame) and ${\bf V}$ is a generalized velocity, |
| 516 |
> |
${\bf V} = |
| 517 |
> |
\left\{{\bf v},{\bf \omega}\right\}$. The right side of |
| 518 |
> |
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a |
| 519 |
> |
system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf |
| 520 |
> |
F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution |
| 521 |
> |
of the system in Newtonian mechanics is typically done in the lab |
| 522 |
> |
frame, it is convenient to handle the dynamics of rigid bodies in |
| 523 |
> |
body-fixed frames. Thus the friction and random forces on each |
| 524 |
> |
substructure are calculated in a body-fixed frame and may converted |
| 525 |
> |
back to the lab frame using that substructure's rotation matrix (${\bf |
| 526 |
> |
Q}$): |
| 527 |
|
\begin{equation} |
| 528 |
< |
F_{r,i}^b (t) = \left( \begin{array}{l} |
| 529 |
< |
f_{r,i}^b (t) \\ |
| 530 |
< |
\tau _{r,i}^b (t) \\ |
| 531 |
< |
\end{array} \right) = - \left( {\begin{array}{*{20}c} |
| 532 |
< |
{\Xi _{R,t} } & {\Xi _{R,c}^T } \\ |
| 533 |
< |
{\Xi _{R,c} } & {\Xi _{R,r} } \\ |
| 534 |
< |
\end{array}} \right)\left( \begin{array}{l} |
| 535 |
< |
v_{R,i}^b (t) \\ |
| 536 |
< |
\omega _i (t) \\ |
| 528 |
> |
{\bf F}_{f,r} = |
| 529 |
> |
\left( \begin{array}{c} |
| 530 |
> |
{\bf f}_{f,r} \\ |
| 531 |
> |
{\bf \tau}_{f,r} |
| 532 |
> |
\end{array} \right) |
| 533 |
> |
= |
| 534 |
> |
\left( \begin{array}{c} |
| 535 |
> |
{\bf Q}^{T} {\bf f}^{~b}_{f,r} \\ |
| 536 |
> |
{\bf Q}^{T} {\bf \tau}^{~b}_{f,r} |
| 537 |
> |
\end{array} \right) |
| 538 |
> |
\end{equation} |
| 539 |
> |
The body-fixed friction force, ${\bf F}_{f}^{~b}$, is proportional to |
| 540 |
> |
the (body-fixed) velocity at the center of resistance |
| 541 |
> |
${\bf v}_{R}^{~b}$ and the angular velocity ${\bf \omega}$ |
| 542 |
> |
\begin{equation} |
| 543 |
> |
{\bf F}_{f}^{~b}(t) = \left( \begin{array}{l} |
| 544 |
> |
{\bf f}_{f}^{~b}(t) \\ |
| 545 |
> |
{\bf \tau}_{f}^{~b}(t) \\ |
| 546 |
> |
\end{array} \right) = - \left( \begin{array}{*{20}c} |
| 547 |
> |
\Xi_{R}^{tt} & \Xi_{R}^{rt} \\ |
| 548 |
> |
\Xi_{R}^{tr} & \Xi_{R}^{rr} \\ |
| 549 |
> |
\end{array} \right)\left( \begin{array}{l} |
| 550 |
> |
{\bf v}_{R}^{~b}(t) \\ |
| 551 |
> |
{\bf \omega}(t) \\ |
| 552 |
|
\end{array} \right), |
| 553 |
|
\end{equation} |
| 554 |
< |
while the random force $F_{r,i}^l$ is a Gaussian stochastic variable |
| 555 |
< |
with zero mean and variance |
| 554 |
> |
while the random force, ${\bf F}_{r}$, is a Gaussian stochastic |
| 555 |
> |
variable with zero mean and variance, |
| 556 |
|
\begin{equation} |
| 557 |
< |
\left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle = |
| 558 |
< |
\left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle = |
| 559 |
< |
2k_B T\Xi _R \delta (t - t'). \label{randomForce} |
| 557 |
> |
\left\langle {{\bf F}_{r}(t) ({\bf F}_{r}(t'))^T } \right\rangle = |
| 558 |
> |
\left\langle {{\bf F}_{r}^{~b} (t) ({\bf F}_{r}^{~b} (t'))^T } \right\rangle = |
| 559 |
> |
2 k_B T \Xi_R \delta(t - t'). \label{eq:randomForce} |
| 560 |
|
\end{equation} |
| 561 |
< |
The equation of motion for $v_i$ can be written as |
| 561 |
> |
$\Xi_R$ is the $6\times6$ resistance tensor at the center of |
| 562 |
> |
resistance. Once this tensor is known for a given rigid body (as |
| 563 |
> |
described in the previous section) obtaining a stochastic vector that |
| 564 |
> |
has the properties in Eq. (\ref{eq:randomForce}) can be done |
| 565 |
> |
efficiently by carrying out a one-time Cholesky decomposition to |
| 566 |
> |
obtain the square root matrix of the resistance tensor, |
| 567 |
> |
\begin{equation} |
| 568 |
> |
\Xi_R = {\bf S} {\bf S}^{T}, |
| 569 |
> |
\label{eq:Cholesky} |
| 570 |
> |
\end{equation} |
| 571 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
| 572 |
> |
vector with the statistics required for the random force can then be |
| 573 |
> |
obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which |
| 574 |
> |
has elements chosen from a Gaussian distribution, such that: |
| 575 |
|
\begin{equation} |
| 576 |
< |
m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) + |
| 577 |
< |
f_{r,i}^l (t) |
| 576 |
> |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
| 577 |
> |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
| 578 |
|
\end{equation} |
| 579 |
< |
Since the frictional force is applied at the center of resistance |
| 580 |
< |
which generally does not coincide with the center of mass, an extra |
| 581 |
< |
torque is exerted at the center of mass. Thus, the net body-fixed |
| 582 |
< |
frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is |
| 583 |
< |
given by |
| 579 |
> |
where $\delta t$ is the timestep in use during the simulation. The |
| 580 |
> |
random force, ${\bf F}_{r}^{~b} = {\bf S} {\bf Z}$, can be shown to have the |
| 581 |
> |
correct properties required by Eq. (\ref{eq:randomForce}). |
| 582 |
> |
|
| 583 |
> |
The equation of motion for the translational velocity of the center of |
| 584 |
> |
mass (${\bf v}$) can be written as |
| 585 |
|
\begin{equation} |
| 586 |
< |
\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b |
| 586 |
> |
m \dot{{\bf v}} (t) = {\bf f}_{s}(t) + {\bf f}_{f}(t) + |
| 587 |
> |
{\bf f}_{r}(t) |
| 588 |
|
\end{equation} |
| 589 |
< |
where $r_{MR}$ is the vector from the center of mass to the center |
| 590 |
< |
of the resistance. Instead of integrating the angular velocity in |
| 591 |
< |
lab-fixed frame, we consider the equation of angular momentum in |
| 592 |
< |
body-fixed frame |
| 589 |
> |
Since the frictional and random forces are applied at the center of |
| 590 |
> |
resistance, which generally does not coincide with the center of mass, |
| 591 |
> |
extra torques are exerted at the center of mass. Thus, the net |
| 592 |
> |
body-fixed torque at the center of mass, $\tau^{~b}(t)$, |
| 593 |
> |
is given by |
| 594 |
|
\begin{equation} |
| 595 |
< |
\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t) |
| 493 |
< |
+ \tau _{r,i}^b(t) |
| 595 |
> |
\tau^{~b} \leftarrow \tau_{s}^{~b} + \tau_{f}^{~b} + \tau_{r}^{~b} + {\bf r}_{MR} \times \left( {\bf f}_{f}^{~b} + {\bf f}_{r}^{~b} \right) |
| 596 |
|
\end{equation} |
| 597 |
< |
Embedding the friction terms into force and torque, one can |
| 598 |
< |
integrate the langevin equations of motion for rigid body of |
| 599 |
< |
arbitrary shape in a velocity-Verlet style 2-part algorithm, where |
| 600 |
< |
$h= \delta t$: |
| 597 |
> |
where ${\bf r}_{MR}$ is the vector from the center of mass to the center of |
| 598 |
> |
resistance. Instead of integrating the angular velocity in lab-fixed |
| 599 |
> |
frame, we consider the equation of motion for the angular momentum |
| 600 |
> |
(${\bf j}$) in the body-fixed frame |
| 601 |
> |
\begin{equation} |
| 602 |
> |
\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 |
| 606 |
> |
body of arbitrary shape in a velocity-Verlet style 2-part algorithm, |
| 607 |
> |
where $h = \delta t$: |
| 608 |
|
|
| 609 |
< |
{\tt moveA:} |
| 609 |
> |
{\tt move A:} |
| 610 |
|
\begin{align*} |
| 611 |
|
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 612 |
|
+ \frac{h}{2} \left( {\bf f}(t) / m \right), \\ |
| 615 |
|
+ h {\bf v}\left(t + h / 2 \right), \\ |
| 616 |
|
% |
| 617 |
|
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 618 |
< |
+ \frac{h}{2} {\bf \tau}^b(t), \\ |
| 618 |
> |
+ \frac{h}{2} {\bf \tau}^{~b}(t), \\ |
| 619 |
|
% |
| 620 |
< |
\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
| 620 |
> |
{\bf Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
| 621 |
|
(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right). |
| 622 |
|
\end{align*} |
| 623 |
< |
In this context, the $\mathrm{rotate}$ function is the reversible |
| 624 |
< |
product of the three body-fixed rotations, |
| 623 |
> |
In this context, $\overleftrightarrow{\mathsf{I}}$ is the diagonal |
| 624 |
> |
moment of inertia tensor, and the $\mathrm{rotate}$ function is the |
| 625 |
> |
reversible product of the three body-fixed rotations, |
| 626 |
|
\begin{equation} |
| 627 |
|
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
| 628 |
|
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y |
| 629 |
|
/ 2) \cdot \mathsf{G}_x(a_x /2), |
| 630 |
|
\end{equation} |
| 631 |
|
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, |
| 632 |
< |
rotates both the rotation matrix ($\mathsf{Q}$) and the body-fixed |
| 632 |
> |
rotates both the rotation matrix ($\mathbf{Q}$) and the body-fixed |
| 633 |
|
angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed |
| 634 |
|
axis $\alpha$, |
| 635 |
|
\begin{equation} |
| 636 |
|
\mathsf{G}_\alpha( \theta ) = \left\{ |
| 637 |
|
\begin{array}{lcl} |
| 638 |
< |
\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
| 638 |
> |
\mathbf{Q}(t) & \leftarrow & \mathbf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
| 639 |
|
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf |
| 640 |
|
j}(0). |
| 641 |
|
\end{array} |
| 657 |
|
\end{equation} |
| 658 |
|
All other rotations follow in a straightforward manner. After the |
| 659 |
|
first part of the propagation, the forces and body-fixed torques are |
| 660 |
< |
calculated at the new positions and orientations |
| 660 |
> |
calculated at the new positions and orientations. The system forces |
| 661 |
> |
and torques are derivatives of the total potential energy function |
| 662 |
> |
($U$) with respect to the rigid body positions (${\bf r}$) and the |
| 663 |
> |
columns of the transposed rotation matrix ${\bf Q}^T = \left({\bf |
| 664 |
> |
u}_x, {\bf u}_y, {\bf u}_z \right)$: |
| 665 |
|
|
| 666 |
< |
{\tt doForces:} |
| 666 |
> |
{\tt Forces:} |
| 667 |
|
\begin{align*} |
| 668 |
< |
{\bf f}(t + h) &\leftarrow |
| 669 |
< |
- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\ |
| 668 |
> |
{\bf f}_{s}(t + h) & \leftarrow |
| 669 |
> |
- \left(\frac{\partial U}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\ |
| 670 |
|
% |
| 671 |
< |
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) |
| 672 |
< |
\times \frac{\partial V}{\partial {\bf u}}, \\ |
| 671 |
> |
{\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h) |
| 672 |
> |
\times \frac{\partial U}{\partial {\bf u}} \\ |
| 673 |
|
% |
| 674 |
< |
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h) |
| 675 |
< |
\cdot {\bf \tau}^s(t + h). |
| 674 |
> |
{\bf v}^{b}_{R}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \left({\bf v}(t+h) + {\bf \omega}(t+h) \times {\bf r}_{MR} \right) \\ |
| 675 |
> |
% |
| 676 |
> |
{\bf f}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tt} \cdot |
| 677 |
> |
{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rt} \cdot {\bf \omega}(t+h) \\ |
| 678 |
> |
% |
| 679 |
> |
{\bf \tau}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tr} \cdot |
| 680 |
> |
{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rr} \cdot {\bf \omega}(t+h) \\ |
| 681 |
> |
% |
| 682 |
> |
Z & \leftarrow {\tt GaussianNormal}(2 k_B T / h, 6) \\ |
| 683 |
> |
{\bf F}_{R,r}^{b}(t+h) & \leftarrow {\bf S} \cdot Z \\ |
| 684 |
> |
% |
| 685 |
> |
{\bf f}(t+h) & \leftarrow {\bf f}_{s}(t+h) + \mathbf{Q}^{T}(t+h) |
| 686 |
> |
\cdot \left({\bf f}_{R,f}^{~b} + {\bf f}_{R,r}^{~b} \right) \\ |
| 687 |
> |
% |
| 688 |
> |
\tau(t+h) & \leftarrow \tau_{s}(t+h) + \mathbf{Q}^{T}(t+h) \cdot \left(\tau_{R,f}^{~b} + \tau_{R,r}^{~b} \right) + {\bf r}_{MR} \times \left({\bf f}_{f}(t+h) + {\bf f}_{r}(t+h) \right) \\ |
| 689 |
> |
\tau^{~b}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \tau(t+h) \\ |
| 690 |
|
\end{align*} |
| 691 |
+ |
Frictional (and random) forces and torques must be computed at the |
| 692 |
+ |
center of resistance, so there are additional steps required to find |
| 693 |
+ |
the body-fixed velocity (${\bf v}_{R}^{~b}$) at this location. Mapping |
| 694 |
+ |
the frictional and random forces at the center of resistance back to |
| 695 |
+ |
the center of mass also introduces an additional term in the torque |
| 696 |
+ |
one obtains at the center of mass. |
| 697 |
+ |
|
| 698 |
|
Once the forces and torques have been obtained at the new time step, |
| 699 |
|
the velocities can be advanced to the same time value. |
| 700 |
|
|
| 701 |
< |
{\tt moveB:} |
| 701 |
> |
{\tt move B:} |
| 702 |
|
\begin{align*} |
| 703 |
|
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2 |
| 704 |
|
\right) |
| 706 |
|
% |
| 707 |
|
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2 |
| 708 |
|
\right) |
| 709 |
< |
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
| 709 |
> |
+ \frac{h}{2} {\bf \tau}^{~b}(t + h) . |
| 710 |
|
\end{align*} |
| 711 |
|
|
| 712 |
< |
\section{Results and Discussion} |
| 712 |
> |
\section{Validating the Method\label{sec:validating}} |
| 713 |
> |
In order to validate our Langevin integrator for arbitrarily-shaped |
| 714 |
> |
rigid bodies, we implemented the algorithm in {\sc |
| 715 |
> |
oopse}\cite{Meineke2005} and compared the results of this algorithm |
| 716 |
> |
with the known |
| 717 |
> |
hydrodynamic limiting behavior for a few model systems, and to |
| 718 |
> |
microcanonical molecular dynamics simulations for some more |
| 719 |
> |
complicated bodies. The model systems and their analytical behavior |
| 720 |
> |
(if known) are summarized below. Parameters for the primary particles |
| 721 |
> |
comprising our model systems are given in table \ref{tab:parameters}, |
| 722 |
> |
and a sketch of the arrangement of these primary particles into the |
| 723 |
> |
model rigid bodies is shown in figure \ref{fig:models}. In table |
| 724 |
> |
\ref{tab:parameters}, $d$ and $l$ are the physical dimensions of |
| 725 |
> |
ellipsoidal (Gay-Berne) particles. For spherical particles, the value |
| 726 |
> |
of the Lennard-Jones $\sigma$ parameter is the particle diameter |
| 727 |
> |
($d$). Gay-Berne ellipsoids have an energy scaling parameter, |
| 728 |
> |
$\epsilon^s$, which describes the well depth for two identical |
| 729 |
> |
ellipsoids in a {\it side-by-side} configuration. Additionally, a |
| 730 |
> |
well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, |
| 731 |
> |
describes the ratio between the well depths in the {\it end-to-end} |
| 732 |
> |
and side-by-side configurations. For spheres, $\epsilon^r \equiv 1$. |
| 733 |
> |
Moments of inertia are also required to describe the motion of primary |
| 734 |
> |
particles with orientational degrees of freedom. |
| 735 |
|
|
| 736 |
< |
The Langevin algorithm described in previous section has been |
| 737 |
< |
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} |
| 736 |
> |
\begin{table*} |
| 737 |
> |
\begin{minipage}{\linewidth} |
| 738 |
|
\begin{center} |
| 739 |
< |
\begin{tabular}{lll} |
| 740 |
< |
\hline |
| 741 |
< |
$\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\ |
| 742 |
< |
\hline |
| 743 |
< |
1 & 300.47 & 10.99 \\ |
| 744 |
< |
0.1 & 301.19 & 11.136 \\ |
| 745 |
< |
0.01 & 303.04 & 11.796 \\ |
| 746 |
< |
\hline |
| 739 |
> |
\caption{Parameters for the primary particles in use by the rigid body |
| 740 |
> |
models in figure \ref{fig:models}.} |
| 741 |
> |
\begin{tabular}{lrcccccccc} |
| 742 |
> |
\hline |
| 743 |
> |
& & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\ |
| 744 |
> |
& & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ (kcal/mol) & $\epsilon^r$ & |
| 745 |
> |
$m$ (amu) & $I_{xx}$ & $I_{yy}$ & $I_{zz}$ \\ \hline |
| 746 |
> |
Sphere & & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
| 747 |
> |
Ellipsoid & & 4.6 & 13.8 & 0.8 & 0.2 & 200 & 2105 & 2105 & 421 \\ |
| 748 |
> |
Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
| 749 |
> |
Banana &(3 identical ellipsoids)& 4.2 & 11.2 & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\ |
| 750 |
> |
Lipid: & Spherical Head & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\ |
| 751 |
> |
& Ellipsoidal Tail & 4.6 & 13.8 & 0.8 & 0.2 & 760 & 45000 & 45000 & 9000 \\ |
| 752 |
> |
Solvent & & 4.7 & $= d$ & 0.8 & 1 & 72.06 & & & \\ |
| 753 |
> |
\hline |
| 754 |
|
\end{tabular} |
| 755 |
+ |
\label{tab:parameters} |
| 756 |
|
\end{center} |
| 757 |
< |
\end{table} |
| 757 |
> |
\end{minipage} |
| 758 |
> |
\end{table*} |
| 759 |
|
|
| 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 |
– |
|
| 760 |
|
\begin{figure} |
| 761 |
|
\centering |
| 762 |
< |
\includegraphics[width=\linewidth]{temperature.pdf} |
| 763 |
< |
\caption[Plot of Temperature Fluctuation Versus Time]{Plot of |
| 764 |
< |
temperature fluctuation versus time.} \label{langevin:temperature} |
| 762 |
> |
\includegraphics[width=3in]{sketch} |
| 763 |
> |
\caption[Sketch of the model systems]{A sketch of the model systems |
| 764 |
> |
used in evaluating the behavior of the rigid body Langevin |
| 765 |
> |
integrator.} \label{fig:models} |
| 766 |
|
\end{figure} |
| 767 |
|
|
| 768 |
< |
\subsection{Langevin Dynamics simulation {\it vs} NVE simulations} |
| 768 |
> |
\subsection{Simulation Methodology} |
| 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 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 (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 |
< |
To validate our langevin dynamics simulation. We performed several NVE |
| 780 |
< |
simulations with explicit solvents for different shaped |
| 781 |
< |
molecules. There are one solute molecule and 1929 solvent molecules in |
| 782 |
< |
NVE simulation. The parameters are shown in table |
| 783 |
< |
\ref{tab:parameters}. The force field between spheres is standard |
| 784 |
< |
Lennard-Jones, and ellipsoids interact with other ellipsoids and |
| 785 |
< |
spheres with generalized Gay-Berne potential. All simulations are |
| 786 |
< |
carried out at 300 K and 1 Atm. The time step is 25 ns, and a |
| 787 |
< |
switching function was applied to all potentials to smoothly turn off |
| 788 |
< |
the interactions between a range of $22$ and $25$ \AA. The switching |
| 789 |
< |
function was the standard (cubic) function, |
| 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{Gay1981} It can be thought of as a |
| 783 |
> |
modification of the Gaussian overlap model originally described by |
| 784 |
> |
Berne and Pechukas.\cite{Berne72} The potential is constructed in the |
| 785 |
> |
familiar form of the Lennard-Jones function using |
| 786 |
> |
orientation-dependent $\sigma$ and $\epsilon$ parameters, |
| 787 |
> |
\begin{equation*} |
| 788 |
> |
V_{ij}({{\bf \hat u}_i}, {{\bf \hat u}_j}, {{\bf \hat |
| 789 |
> |
r}_{ij}}) = 4\epsilon ({{\bf \hat u}_i}, {{\bf \hat u}_j}, |
| 790 |
> |
{{\bf \hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u |
| 791 |
> |
}_i}, |
| 792 |
> |
{{\bf \hat u}_j}, {{\bf \hat r}_{ij}})+\sigma_0}\right)^{12} |
| 793 |
> |
-\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u}_i}, {{\bf \hat u}_j}, |
| 794 |
> |
{{\bf \hat r}_{ij}})+\sigma_0}\right)^6\right] |
| 795 |
> |
\label{eq:gb} |
| 796 |
> |
\end{equation*} |
| 797 |
> |
|
| 798 |
> |
The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf |
| 799 |
> |
\hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf |
| 800 |
> |
\hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters |
| 801 |
> |
are dependent on the relative orientations of the two ellipsoids (${\bf |
| 802 |
> |
\hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the |
| 803 |
> |
inter-ellipsoid separation (${\bf \hat{r}}_{ij}$). The shape and |
| 804 |
> |
attractiveness of each ellipsoid is governed by a relatively small set |
| 805 |
> |
of parameters: $l$ and $d$ describe the length and width of each |
| 806 |
> |
uniaxial ellipsoid, while $\epsilon^s$, which describes the well depth |
| 807 |
> |
for two identical ellipsoids in a {\it side-by-side} configuration. |
| 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,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} |
| 815 |
> |
|
| 816 |
> |
For the interaction between nonequivalent uniaxial ellipsoids (or |
| 817 |
> |
between spheres and ellipsoids), the spheres are treated as ellipsoids |
| 818 |
> |
with an aspect ratio of 1 ($d = l$) and with an well depth ratio |
| 819 |
> |
($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$). The form of the |
| 820 |
> |
Gay-Berne potential we are using was generalized by Cleaver {\it et |
| 821 |
> |
al.} and is appropriate for dissimilar uniaxial |
| 822 |
> |
ellipsoids.\cite{Cleaver96} |
| 823 |
> |
|
| 824 |
> |
A switching function was applied to all potentials to smoothly turn |
| 825 |
> |
off the interactions between a range of $22$ and $25$ \AA. The |
| 826 |
> |
switching function was the standard (cubic) function, |
| 827 |
|
\begin{equation} |
| 828 |
|
s(r) = |
| 829 |
|
\begin{cases} |
| 835 |
|
\end{cases} |
| 836 |
|
\label{eq:switchingFunc} |
| 837 |
|
\end{equation} |
| 838 |
< |
We have computed translational diffusion constants for lipid molecules |
| 839 |
< |
from the mean-square displacement, |
| 838 |
> |
|
| 839 |
> |
To measure shear viscosities from our microcanonical simulations, we |
| 840 |
> |
used the Einstein form of the pressure correlation function,\cite{hess:209} |
| 841 |
|
\begin{equation} |
| 842 |
< |
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle, |
| 842 |
> |
\eta = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left( |
| 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 |
< |
of the solute molecules. Translational diffusion constants for the |
| 847 |
< |
different shaped molecules are shown in table |
| 848 |
< |
\ref{tab:translation}. We have also computed orientational correlation |
| 849 |
< |
times for different shaped molecules from fits of the second-order |
| 850 |
< |
Legendre polynomial correlation function, |
| 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 |
| 851 |
> |
microcanonical simulations and with a solvent viscosity taken from |
| 852 |
> |
Eq. (\ref{eq:shear}) applied to these simulations. We used 1024 |
| 853 |
> |
independent solute simulations to obtain statistics on our Langevin |
| 854 |
> |
integrator. |
| 855 |
> |
|
| 856 |
> |
\subsection{Analysis} |
| 857 |
> |
|
| 858 |
> |
The quantities of interest when comparing the Langevin integrator to |
| 859 |
> |
analytic hydrodynamic equations and to molecular dynamics simulations |
| 860 |
> |
are typically translational diffusion constants and orientational |
| 861 |
> |
relaxation times. Translational diffusion constants for point |
| 862 |
> |
particles are computed easily from the long-time slope of the |
| 863 |
> |
mean-square displacement, |
| 864 |
|
\begin{equation} |
| 865 |
< |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
| 675 |
< |
\mu}_{i}(0) \right) |
| 865 |
> |
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \left\langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \right\rangle, |
| 866 |
|
\end{equation} |
| 867 |
< |
the results are shown in table \ref{tab:rotation}. We used einstein |
| 868 |
< |
format of the pressure correlation function, |
| 867 |
> |
of the solute molecules. For models in which the translational |
| 868 |
> |
diffusion tensor (${\bf D}_{tt}$) has non-degenerate eigenvalues |
| 869 |
> |
(i.e. any non-spherically-symmetric rigid body), it is possible to |
| 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 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 |
| 880 |
> |
{\it around} a particular body-fixed axis and {\it not} the diffusion |
| 881 |
> |
of a vector pointing along the axis. However, these eigenvalues can |
| 882 |
> |
be combined to find 5 characteristic rotational relaxation |
| 883 |
> |
times,\cite{PhysRev.119.53,Berne90} |
| 884 |
> |
\begin{eqnarray} |
| 885 |
> |
1 / \tau_1 & = & 6 D_r + 2 \Delta \\ |
| 886 |
> |
1 / \tau_2 & = & 6 D_r - 2 \Delta \\ |
| 887 |
> |
1 / \tau_3 & = & 3 (D_r + D_1) \\ |
| 888 |
> |
1 / \tau_4 & = & 3 (D_r + D_2) \\ |
| 889 |
> |
1 / \tau_5 & = & 3 (D_r + D_3) |
| 890 |
> |
\end{eqnarray} |
| 891 |
> |
where |
| 892 |
|
\begin{equation} |
| 893 |
< |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
| 681 |
< |
\mu}_{i}(0) \right) |
| 893 |
> |
D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right) |
| 894 |
|
\end{equation} |
| 895 |
< |
to estimate the viscosity of the systems from NVE simulations. The |
| 684 |
< |
viscosity can also be calculated by Green-Kubo pressure correlaton |
| 685 |
< |
function, |
| 895 |
> |
and |
| 896 |
|
\begin{equation} |
| 897 |
< |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
| 688 |
< |
\mu}_{i}(0) \right) |
| 897 |
> |
\Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2} |
| 898 |
|
\end{equation} |
| 899 |
< |
However, this method converges slowly, and the statistics are not good |
| 900 |
< |
enough to give us a very accurate value. The langevin dynamics |
| 901 |
< |
simulations for different shaped molecules are performed at the same |
| 902 |
< |
conditions as the NVE simulations with viscosity estimated from NVE |
| 903 |
< |
simulations. To get better statistics, 1024 non-interacting solute |
| 904 |
< |
molecules are put into one simulation box for each langevin |
| 905 |
< |
simulation, this is equal to 1024 simulations for single solute |
| 906 |
< |
systems. The diffusion constants and rotation relaxation times for |
| 907 |
< |
different shaped molecules are shown in table \ref{tab:translation} |
| 908 |
< |
and \ref{tab:rotation} to compare to the results calculated from NVE |
| 909 |
< |
simulations. The theoretical values for sphere is calculated from the |
| 910 |
< |
Stokes-Einstein law, the theoretical values for ellipsoid is |
| 911 |
< |
calculated from Perrin's fomula, the theoretical values for dumbbell |
| 912 |
< |
is calculated from StinXX and Davis theory. The exact method is |
| 913 |
< |
applied to the langevin dynamics simulations for sphere and ellipsoid, |
| 914 |
< |
the bead model is applied to the simulation for dumbbell molecule, and |
| 915 |
< |
the rough shell model is applied to ellipsoid, dumbbell, banana and |
| 916 |
< |
lipid molecules. The results from all the langevin dynamics |
| 917 |
< |
simulations, including exact, bead model and rough shell, match the |
| 918 |
< |
theoretical values perfectly for all different shaped molecules. This |
| 919 |
< |
indicates that our simulation package for langevin dynamics is working |
| 920 |
< |
well. The approxiate methods ( bead model and rough shell model) are |
| 712 |
< |
accurate enough for the current simulations. The goal of the langevin |
| 713 |
< |
dynamics theory is to replace the explicit solvents by the friction |
| 714 |
< |
forces. We compared the dynamic properties of different shaped |
| 715 |
< |
molecules in langevin dynamics simulations with that in NVE |
| 716 |
< |
simulations. The results are reasonable close. Overall, the |
| 717 |
< |
translational diffusion constants calculated from langevin dynamics |
| 718 |
< |
simulations are very close to the values from the NVE simulation. For |
| 719 |
< |
sphere and lipid molecules, the diffusion constants are a little bit |
| 720 |
< |
off from the NVE simulation results. One possible reason is that the |
| 721 |
< |
calculation of the viscosity is very difficult to be accurate. Another |
| 722 |
< |
possible reason is that although we save very frequently during the |
| 723 |
< |
NVE simulations and run pretty long time simulations, there is only |
| 724 |
< |
one solute molecule in the system which makes the calculation for the |
| 725 |
< |
diffusion constant difficult. The sphere molecule behaves as a free |
| 726 |
< |
rotor in the solvent, so there is no rotation relaxation time |
| 727 |
< |
calculated from NVE simulations. The rotation relaxation time is not |
| 728 |
< |
very close to the NVE simulations results. The banana and lipid |
| 729 |
< |
molecules match the NVE simulations results pretty well. The mismatch |
| 730 |
< |
between langevin dynamics and NVE simulation for ellipsoid is possibly |
| 731 |
< |
caused by the slip boundary condition. For dumbbell, the mismatch is |
| 732 |
< |
caused by the size of the solvent molecule is pretty large compared to |
| 733 |
< |
dumbbell molecule in NVE simulations. |
| 899 |
> |
Each of these characteristic times can be used to predict the decay of |
| 900 |
> |
part of the rotational correlation function when $\ell = 2$, |
| 901 |
> |
\begin{equation} |
| 902 |
> |
C_2(t) = \frac{a^2}{N^2} e^{-t/\tau_1} + \frac{b^2}{N^2} e^{-t/\tau_2}. |
| 903 |
> |
\end{equation} |
| 904 |
> |
This is the same as the $F^2_{0,0}(t)$ correlation function that |
| 905 |
> |
appears in Ref. \citen{Berne90}. The amplitudes of the two decay |
| 906 |
> |
terms are expressed in terms of three dimensionless functions of the |
| 907 |
> |
eigenvalues: $a = \sqrt{3} (D_1 - D_2)$, $b = (2D_3 - D_1 - D_2 + |
| 908 |
> |
2\Delta)$, and $N = 2 \sqrt{\Delta b}$. Similar expressions can be |
| 909 |
> |
obtained for other angular momentum correlation |
| 910 |
> |
functions.\cite{PhysRev.119.53,Berne90} In all of the model systems we |
| 911 |
> |
studied, only one of the amplitudes of the two decay terms was |
| 912 |
> |
non-zero, so it was possible to derive a single relaxation time for |
| 913 |
> |
each of the hydrodynamic tensors. In many cases, these characteristic |
| 914 |
> |
times are averaged and reported in the literature as a single relaxation |
| 915 |
> |
time,\cite{Garcia-de-la-Torre:1997qy} |
| 916 |
> |
\begin{equation} |
| 917 |
> |
1 / \tau_0 = \frac{1}{5} \sum_{i=1}^5 \tau_{i}^{-1}, |
| 918 |
> |
\end{equation} |
| 919 |
> |
although for the cases reported here, this averaging is not necessary |
| 920 |
> |
and only one of the five relaxation times is relevant. |
| 921 |
|
|
| 922 |
< |
According to our simulations, the langevin dynamics is a reliable |
| 923 |
< |
theory to apply to replace the explicit solvents, especially for the |
| 924 |
< |
translation properties. For large molecules, the rotation properties |
| 925 |
< |
are also mimiced reasonablly well. |
| 922 |
> |
To test the Langevin integrator's behavior for rotational relaxation, |
| 923 |
> |
we have compared the analytical orientational relaxation times (if |
| 924 |
> |
they are known) with the general result from the diffusion tensor and |
| 925 |
> |
with the results from both the explicitly solvated molecular dynamics |
| 926 |
> |
and Langevin simulations. Relaxation times from simulations (both |
| 927 |
> |
microcanonical and Langevin), were computed using Legendre polynomial |
| 928 |
> |
correlation functions for a unit vector (${\bf u}$) fixed along one or |
| 929 |
> |
more of the body-fixed axes of the model. |
| 930 |
> |
\begin{equation} |
| 931 |
> |
C_{\ell}(t) = \left\langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf |
| 932 |
> |
u}_{i}(0) \right) \right\rangle |
| 933 |
> |
\end{equation} |
| 934 |
> |
For simulations in the high-friction limit, orientational correlation |
| 935 |
> |
times can then be obtained from exponential fits of this function, or by |
| 936 |
> |
integrating, |
| 937 |
> |
\begin{equation} |
| 938 |
> |
\tau = \ell (\ell + 1) \int_0^{\infty} C_{\ell}(t) dt. |
| 939 |
> |
\end{equation} |
| 940 |
> |
In lower-friction solvents, the Legendre correlation functions often |
| 941 |
> |
exhibit non-exponential decay, and may not be characterized by a |
| 942 |
> |
single decay constant. |
| 943 |
|
|
| 944 |
< |
\begin{table*} |
| 945 |
< |
\begin{minipage}{\linewidth} |
| 946 |
< |
\begin{center} |
| 947 |
< |
\caption{} |
| 948 |
< |
\begin{tabular}{llccccccc} |
| 949 |
< |
\hline |
| 950 |
< |
& & Sphere & Ellipsoid & Dumbbell(2 spheres) & Banana(3 ellpsoids) & |
| 951 |
< |
Lipid(head) & lipid(tail) & Solvent \\ |
| 952 |
< |
\hline |
| 953 |
< |
$d$ (\AA) & & 6.5 & 4.6 & 6.5 & 4.2 & 6.5 & 4.6 & 4.7 \\ |
| 954 |
< |
$l$ (\AA) & & $= d$ & 13.8 & $=d$ & 11.2 & $=d$ & 13.8 & 4.7 \\ |
| 955 |
< |
$\epsilon^s$ (kcal/mol) & & 0.8 & 0.8 & 0.8 & 0.8 & 0.185 & 0.8 & 0.8 \\ |
| 956 |
< |
$\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 & 1 & 0.2 & 1 & 0.2 & 1 \\ |
| 957 |
< |
$m$ (amu) & & 190 & 200 & 190 & 240 & 196 & 760 & 72.06 \\ |
| 958 |
< |
%$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\ |
| 959 |
< |
%\multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\ |
| 960 |
< |
%\multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\ |
| 961 |
< |
%\multicolumn{2}c{$I_{zz}$} & 0 & 9000 & N/A \\ |
| 962 |
< |
%$\mu$ (Debye) & & varied & 0 & 0 \\ |
| 963 |
< |
\end{tabular} |
| 964 |
< |
\label{tab:parameters} |
| 965 |
< |
\end{center} |
| 966 |
< |
\end{minipage} |
| 967 |
< |
\end{table*} |
| 944 |
> |
In table \ref{tab:rotation} we show the characteristic rotational |
| 945 |
> |
relaxation times (based on the diffusion tensor) for each of the model |
| 946 |
> |
systems compared with the values obtained via microcanonical and Langevin |
| 947 |
> |
simulations. |
| 948 |
> |
|
| 949 |
> |
\subsection{Spherical particles} |
| 950 |
> |
Our model system for spherical particles was a Lennard-Jones sphere of |
| 951 |
> |
diameter ($\sigma$) 6.5 \AA\ in a sea of smaller spheres ($\sigma$ = |
| 952 |
> |
4.7 \AA). The well depth ($\epsilon$) for both particles was set to |
| 953 |
> |
an arbitrary value of 0.8 kcal/mol. |
| 954 |
> |
|
| 955 |
> |
The Stokes-Einstein behavior of large spherical particles in |
| 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 |
| 967 |
> |
microcanonical molecular dynamics simulations lies between the slip |
| 968 |
> |
and stick boundary condition results obtained via Stokes-Einstein |
| 969 |
> |
behavior. Since the Langevin integrator assumes Stokes-Einstein stick |
| 970 |
> |
boundary conditions in calculating the drag and random forces for |
| 971 |
> |
spherical particles, our Langevin routine obtains nearly quantitative |
| 972 |
> |
agreement with the hydrodynamic results for spherical particles. One |
| 973 |
> |
avenue for improvement of the method would be to compute elements of |
| 974 |
> |
$\Xi_{tt}$ assuming behavior intermediate between the two boundary |
| 975 |
> |
conditions. |
| 976 |
> |
|
| 977 |
> |
In the explicit solvent simulations, both our solute and solvent |
| 978 |
> |
particles were structureless, exerting no torques upon each other. |
| 979 |
> |
Therefore, there are not rotational correlation times available for |
| 980 |
> |
this model system. |
| 981 |
> |
|
| 982 |
> |
\subsection{Ellipsoids} |
| 983 |
> |
For uniaxial ellipsoids ($a > b = c$), Perrin's formulae for both |
| 984 |
> |
translational and rotational diffusion of each of the body-fixed axes |
| 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(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 - s^2) |
| 994 |
> |
G(s) - 1}{1 - s^4} \right\}. |
| 995 |
> |
\label{ThetaPerrin} |
| 996 |
> |
\end{equation} |
| 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(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 |
| 1004 |
> |
to use for molecular-scale ellipsoids in a sea of similarly-sized |
| 1005 |
> |
solvent particles. Ravichandran and Bagchi found that {\it slip} |
| 1006 |
> |
boundary conditions most closely resembled the simulation |
| 1007 |
> |
results,\cite{Ravichandran:1999fk} in agreement with earlier work of |
| 1008 |
> |
Tang and Evans.\cite{TANG:1993lr} |
| 1009 |
|
|
| 1010 |
+ |
Even though there are analytic resistance tensors for ellipsoids, we |
| 1011 |
+ |
constructed a rough-shell model using 2135 beads (each with a diameter |
| 1012 |
+ |
of 0.25 \AA) to approximate the shape of the model ellipsoid. We |
| 1013 |
+ |
compared the Langevin dynamics from both the simple ellipsoidal |
| 1014 |
+ |
resistance tensor and the rough shell approximation with |
| 1015 |
+ |
microcanonical simulations and the predictions of Perrin. As in the |
| 1016 |
+ |
case of our spherical model system, the Langevin integrator reproduces |
| 1017 |
+ |
almost exactly the behavior of the Perrin formulae (which is |
| 1018 |
+ |
unsurprising given that the Perrin formulae were used to derive the |
| 1019 |
+ |
drag and random forces applied to the ellipsoid). We obtain |
| 1020 |
+ |
translational diffusion constants and rotational correlation times |
| 1021 |
+ |
that are within a few percent of the analytic values for both the |
| 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 |
| 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 |
| 1040 |
+ |
expressions for the hydrodynamic tensor are available is the |
| 1041 |
+ |
two-sphere dumbbell model. This model consists of two non-overlapping |
| 1042 |
+ |
spheres held by a rigid bond connecting their centers. There are |
| 1043 |
+ |
competing expressions for the 6x6 resistance tensor for this |
| 1044 |
+ |
model. The second order expression introduced by Rotne and |
| 1045 |
+ |
Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la Torre and |
| 1046 |
+ |
Bloomfield,\cite{Torre1977} is given above as |
| 1047 |
+ |
Eq. (\ref{introEquation:RPTensorNonOverlapped}). In our case, we use |
| 1048 |
+ |
a model dumbbell in which the two spheres are identical Lennard-Jones |
| 1049 |
+ |
particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at |
| 1050 |
+ |
a distance of 6.532 \AA. |
| 1051 |
+ |
|
| 1052 |
+ |
The theoretical values for the translational diffusion constant of the |
| 1053 |
+ |
dumbbell are calculated from the work of Stimson and Jeffery, who |
| 1054 |
+ |
studied the motion of this system in a flow parallel to the |
| 1055 |
+ |
inter-sphere axis,\cite{Stimson:1926qy} and Davis, who studied the |
| 1056 |
+ |
motion in a flow {\it perpendicular} to the inter-sphere |
| 1057 |
+ |
axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it |
| 1058 |
+ |
orientational} correlation times for this model system (other than |
| 1059 |
+ |
those derived from the 6 x 6 tensor mentioned above). |
| 1060 |
+ |
|
| 1061 |
+ |
The bead model for this model system comprises the two large spheres |
| 1062 |
+ |
by themselves, while the rough shell approximation used 3368 separate |
| 1063 |
+ |
beads (each with a diameter of 0.25 \AA) to approximate the shape of |
| 1064 |
+ |
the rigid body. The hydrodynamics tensors computed from both the bead |
| 1065 |
+ |
and rough shell models are remarkably similar. Computing the initial |
| 1066 |
+ |
hydrodynamic tensor for a rough shell model can be quite expensive (in |
| 1067 |
+ |
this case it requires inverting a 10104 x 10104 matrix), while the |
| 1068 |
+ |
bead model is typically easy to compute (in this case requiring |
| 1069 |
+ |
inversion of a 6 x 6 matrix). |
| 1070 |
+ |
|
| 1071 |
+ |
\begin{figure} |
| 1072 |
+ |
\centering |
| 1073 |
+ |
\includegraphics[width=2in]{RoughShell} |
| 1074 |
+ |
\caption[Model rigid bodies and their rough shell approximations]{The |
| 1075 |
+ |
model rigid bodies (left column) used to test this algorithm and their |
| 1076 |
+ |
rough-shell approximations (right-column) that were used to compute |
| 1077 |
+ |
the hydrodynamic tensors. The top two models (ellipsoid and dumbbell) |
| 1078 |
+ |
have analytic solutions and were used to test the rough shell |
| 1079 |
+ |
approximation. The lower two models (banana and lipid) were compared |
| 1080 |
+ |
with explicitly-solvated molecular dynamics simulations. } |
| 1081 |
+ |
\label{fig:roughShell} |
| 1082 |
+ |
\end{figure} |
| 1083 |
+ |
|
| 1084 |
+ |
|
| 1085 |
+ |
Once the hydrodynamic tensor has been computed, there is no additional |
| 1086 |
+ |
penalty for carrying out a Langevin simulation with either of the two |
| 1087 |
+ |
different hydrodynamics models. Our naive expectation is that since |
| 1088 |
+ |
the rigid body's surface is roughened under the various shell models, |
| 1089 |
+ |
the diffusion constants will be even farther from the ``slip'' |
| 1090 |
+ |
boundary conditions than observed for the bead model (which uses a |
| 1091 |
+ |
Stokes-Einstein model to arrive at the hydrodynamic tensor). For the |
| 1092 |
+ |
dumbbell, this prediction is correct although all of the Langevin |
| 1093 |
+ |
diffusion constants are within 6\% of the diffusion constant predicted |
| 1094 |
+ |
from the fully solvated system. |
| 1095 |
+ |
|
| 1096 |
+ |
For rotational motion, Langevin integration (and the hydrodynamic tensor) |
| 1097 |
+ |
yields rotational correlation times that are substantially shorter than those |
| 1098 |
+ |
obtained from explicitly-solvated simulations. It is likely that this is due |
| 1099 |
+ |
to the large size of the explicit solvent spheres, a feature that prevents |
| 1100 |
+ |
the solvent from coming in contact with a substantial fraction of the surface |
| 1101 |
+ |
area of the dumbbell. Therefore, the explicit solvent only provides drag |
| 1102 |
+ |
over a substantially reduced surface area of this model, while the |
| 1103 |
+ |
hydrodynamic theories utilize the entire surface area for estimating |
| 1104 |
+ |
rotational diffusion. A sketch of the free volume available in the explicit |
| 1105 |
+ |
solvent simulations is shown in figure \ref{fig:explanation}. |
| 1106 |
+ |
|
| 1107 |
+ |
|
| 1108 |
+ |
\begin{figure} |
| 1109 |
+ |
\centering |
| 1110 |
+ |
\includegraphics[width=6in]{explanation} |
| 1111 |
+ |
\caption[Explanations of the differences between orientational |
| 1112 |
+ |
correlation times for explicitly-solvated models and hydrodynamics |
| 1113 |
+ |
predictions]{Explanations of the differences between orientational |
| 1114 |
+ |
correlation times for explicitly-solvated models and hydrodynamic |
| 1115 |
+ |
predictions. For the ellipsoids (upper figures), rotation of the |
| 1116 |
+ |
principal axis can involve correlated collisions at both sides of the |
| 1117 |
+ |
solute. In the rigid dumbbell model (lower figures), the large size |
| 1118 |
+ |
of the explicit solvent spheres prevents them from coming in contact |
| 1119 |
+ |
with a substantial fraction of the surface area of the dumbbell. |
| 1120 |
+ |
Therefore, the explicit solvent only provides drag over a |
| 1121 |
+ |
substantially reduced surface area of this model, where the |
| 1122 |
+ |
hydrodynamic theories utilize the entire surface area for estimating |
| 1123 |
+ |
rotational diffusion. |
| 1124 |
+ |
} \label{fig:explanation} |
| 1125 |
+ |
\end{figure} |
| 1126 |
+ |
|
| 1127 |
+ |
\subsection{Composite banana-shaped molecules} |
| 1128 |
+ |
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids have |
| 1129 |
+ |
been used by Orlandi {\it et al.} to observe mesophases in |
| 1130 |
+ |
coarse-grained models for bent-core liquid crystalline |
| 1131 |
+ |
molecules.\cite{Orlandi:2006fk} We have used the same overlapping |
| 1132 |
+ |
ellipsoids as a way to test the behavior of our algorithm for a |
| 1133 |
+ |
structure of some interest to the materials science community, |
| 1134 |
+ |
although since we are interested in capturing only the hydrodynamic |
| 1135 |
+ |
behavior of this model, we have left out the dipolar interactions of |
| 1136 |
+ |
the original Orlandi model. |
| 1137 |
+ |
|
| 1138 |
+ |
A reference system composed of a single banana rigid body embedded in |
| 1139 |
+ |
a sea of 1929 solvent particles was created and run under standard |
| 1140 |
+ |
(microcanonical) molecular dynamics. The resulting viscosity of this |
| 1141 |
+ |
mixture was 0.298 centipoise (as estimated using |
| 1142 |
+ |
Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties of |
| 1143 |
+ |
the banana rigid body model, we created a rough shell (see |
| 1144 |
+ |
Fig.~\ref{fig:roughShell}), in which the banana is represented as a |
| 1145 |
+ |
``shell'' made of 3321 identical beads (0.25 \AA\ in diameter) |
| 1146 |
+ |
distributed on the surface. Applying the procedure described in |
| 1147 |
+ |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
| 1148 |
+ |
identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 |
| 1149 |
+ |
\AA). |
| 1150 |
+ |
|
| 1151 |
+ |
The Langevin rigid-body integrator (and the hydrodynamic diffusion |
| 1152 |
+ |
tensor) are essentially quantitative for translational diffusion of |
| 1153 |
+ |
this model. Orientational correlation times under the Langevin |
| 1154 |
+ |
rigid-body integrator are within 11\% of the values obtained from |
| 1155 |
+ |
explicit solvent, but these models also exhibit some solvent |
| 1156 |
+ |
inaccessible surface area in the explicitly-solvated case. |
| 1157 |
+ |
|
| 1158 |
+ |
\subsection{Composite sphero-ellipsoids} |
| 1159 |
+ |
|
| 1160 |
+ |
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
| 1161 |
+ |
used recently as models for lipid |
| 1162 |
+ |
molecules.\cite{SunX._jp0762020,Ayton01} A reference system composed of |
| 1163 |
+ |
a single lipid rigid body embedded in a sea of 1929 solvent particles |
| 1164 |
+ |
was created and run under a microcanonical ensemble. The resulting |
| 1165 |
+ |
viscosity of this mixture was 0.349 centipoise (as estimated using |
| 1166 |
+ |
Eq. (\ref{eq:shear})). To calculate the hydrodynamic properties of |
| 1167 |
+ |
the lipid rigid body model, we created a rough shell (see |
| 1168 |
+ |
Fig.~\ref{fig:roughShell}), in which the lipid is represented as a |
| 1169 |
+ |
``shell'' made of 3550 identical beads (0.25 \AA\ in diameter) |
| 1170 |
+ |
distributed on the surface. Applying the procedure described in |
| 1171 |
+ |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
| 1172 |
+ |
identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46 |
| 1173 |
+ |
\AA). |
| 1174 |
+ |
|
| 1175 |
+ |
The translational diffusion constants and rotational correlation times |
| 1176 |
+ |
obtained using the Langevin rigid-body integrator (and the |
| 1177 |
+ |
hydrodynamic tensor) are essentially quantitative when compared with |
| 1178 |
+ |
the explicit solvent simulations for this model system. |
| 1179 |
+ |
|
| 1180 |
+ |
\subsection{Summary 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 |
| 1196 |
+ |
\includegraphics[width=\linewidth]{graph} |
| 1197 |
+ |
\caption[Mean squared displacements and orientational |
| 1198 |
+ |
correlation functions for each of the model rigid bodies.]{The |
| 1199 |
+ |
mean-squared displacements ($\langle r^2(t) \rangle$) and |
| 1200 |
+ |
orientational correlation functions ($C_2(t)$) for each of the model |
| 1201 |
+ |
rigid bodies studied. The circles are the results for microcanonical |
| 1202 |
+ |
simulations with explicit solvent molecules, while the other data sets |
| 1203 |
+ |
are results for Langevin dynamics using the different hydrodynamic |
| 1204 |
+ |
tensor approximations. The Perrin model for the ellipsoids is |
| 1205 |
+ |
considered the ``exact'' hydrodynamic behavior (this can also be said |
| 1206 |
+ |
for the translational motion of the dumbbell operating under the bead |
| 1207 |
+ |
model). In most cases, the various hydrodynamics models reproduce |
| 1208 |
+ |
each other quantitatively.} |
| 1209 |
+ |
\label{fig:results} |
| 1210 |
+ |
\end{figure} |
| 1211 |
+ |
|
| 1212 |
|
\begin{table*} |
| 1213 |
|
\begin{minipage}{\linewidth} |
| 1214 |
|
\begin{center} |
| 1215 |
< |
\caption{} |
| 1216 |
< |
\begin{tabular}{lccccc} |
| 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 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 (= |
| 1223 |
> |
$10^{-4}$ \AA$^2$ / fs). } |
| 1224 |
> |
\begin{tabular}{lccccccc} |
| 1225 |
|
\hline |
| 1226 |
< |
& & & & &Translation \\ |
| 1227 |
< |
\hline |
| 1228 |
< |
& NVE & & Theoretical & Langevin & \\ |
| 774 |
< |
\hline |
| 775 |
< |
& $\eta$ & D & D & method & D \\ |
| 1226 |
> |
& \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\ |
| 1227 |
> |
\cline{2-3} \cline{5-7} |
| 1228 |
> |
model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\ |
| 1229 |
|
\hline |
| 1230 |
< |
sphere & 3.480159e-03 & 1.643135e-04 & 1.942779e-04 & exact & 1.982283e-04 \\ |
| 1231 |
< |
ellipsoid & 2.551262e-03 & 2.437492e-04 & 2.335756e-04 & exact & 2.374905e-04 \\ |
| 1232 |
< |
& 2.551262e-03 & 2.437492e-04 & 2.335756e-04 & rough shell & 2.284088e-04 \\ |
| 1233 |
< |
dumbell & 2.41276e-03 & 2.129432e-04 & 2.090239e-04 & bead model & 2.148098e-04 \\ |
| 1234 |
< |
& 2.41276e-03 & 2.129432e-04 & 2.090239e-04 & rough shell & 2.013219e-04 \\ |
| 1235 |
< |
banana & 2.9846e-03 & 1.527819e-04 & & rough shell & 1.54807e-04 \\ |
| 1236 |
< |
lipid & 3.488661e-03 & 0.9562979e-04 & & rough shell & 1.320987e-04 \\ |
| 1230 |
> |
sphere & 0.279 & 3.06 & & 2.42 & exact & 2.42 & 2.33 \\ |
| 1231 |
> |
ellipsoid & 0.255 & 2.44 & & 2.34 & exact & 2.34 & 2.37 \\ |
| 1232 |
> |
& 0.255 & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\ |
| 1233 |
> |
dumbbell & 0.308 & 2.06 & & 1.64 & bead model & 1.65 & 1.62 \\ |
| 1234 |
> |
& 0.308 & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\ |
| 1235 |
> |
banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\ |
| 1236 |
> |
lipid & 0.349 & 1.41 & & & rough shell & 1.33 & 1.32 \\ |
| 1237 |
|
\end{tabular} |
| 1238 |
|
\label{tab:translation} |
| 1239 |
|
\end{center} |
| 1243 |
|
\begin{table*} |
| 1244 |
|
\begin{minipage}{\linewidth} |
| 1245 |
|
\begin{center} |
| 1246 |
< |
\caption{} |
| 1247 |
< |
\begin{tabular}{lccccc} |
| 1246 |
> |
\caption{Orientational relaxation times ($\tau$) for the model systems using |
| 1247 |
> |
microcanonical simulation (with explicit solvent), theoretical |
| 1248 |
> |
predictions, and Langevin simulations (with implicit solvent). All |
| 1249 |
> |
relaxation times are for the rotational correlation function with |
| 1250 |
> |
$\ell = 2$ and are reported in units of ps. The ellipsoidal model has |
| 1251 |
> |
an exact solution for the orientational correlation time due to |
| 1252 |
> |
Perrin, but the other model systems have no known analytic solution.} |
| 1253 |
> |
\begin{tabular}{lccccccc} |
| 1254 |
|
\hline |
| 1255 |
< |
& & & & &Rotation \\ |
| 1256 |
< |
\hline |
| 1257 |
< |
& NVE & & Theoretical & Langevin & \\ |
| 799 |
< |
\hline |
| 800 |
< |
& $\eta$ & $\tau_0$ & $\tau_0$ & method & $\tau_0$ \\ |
| 1255 |
> |
& \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\ |
| 1256 |
> |
\cline{2-3} \cline{5-7} |
| 1257 |
> |
model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic & simulation \\ |
| 1258 |
|
\hline |
| 1259 |
< |
sphere & 3.480159e-03 & & 1.208178e+04 & exact & 1.20628e+04 \\ |
| 1260 |
< |
ellipsoid & 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & exact & 2.21507e+04 \\ |
| 1261 |
< |
& 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & rough shell & 2.21714e+04 \\ |
| 1262 |
< |
dumbell & 2.41276e-03 & 1.42974e+04 & & bead model & 7.12435e+04 \\ |
| 1263 |
< |
& 2.41276e-03 & 1.42974e+04 & & rough shell & 7.04765e+04 \\ |
| 1264 |
< |
banana & 2.9846e-03 & 6.38323e+04 & & rough shell & 7.0945e+04 \\ |
| 1265 |
< |
lipid & 3.488661e-03 & 7.79595e+04 & & rough shell & 7.78886e+04 \\ |
| 1259 |
> |
sphere & 0.279 & & & 9.69 & exact & 9.69 & 9.64 \\ |
| 1260 |
> |
ellipsoid & 0.255 & 46.7 & & 22.0 & exact & 22.0 & 22.2 \\ |
| 1261 |
> |
& 0.255 & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\ |
| 1262 |
> |
dumbbell & 0.308 & 14.1 & & & bead model & 50.0 & 50.1 \\ |
| 1263 |
> |
& 0.308 & 14.1 & & & rough shell & 41.5 & 41.3 \\ |
| 1264 |
> |
banana & 0.298 & 63.8 & & & rough shell & 70.9 & 70.9 \\ |
| 1265 |
> |
lipid & 0.349 & 78.0 & & & rough shell & 76.9 & 77.9 \\ |
| 1266 |
> |
\hline |
| 1267 |
|
\end{tabular} |
| 1268 |
|
\label{tab:rotation} |
| 1269 |
|
\end{center} |
| 1270 |
|
\end{minipage} |
| 1271 |
|
\end{table*} |
| 1272 |
|
|
| 1273 |
< |
Langevin dynamics simulations are applied to study the formation of |
| 816 |
< |
the ripple phase of lipid membranes. The initial configuration is |
| 817 |
< |
taken from our molecular dynamics studies on lipid bilayers with |
| 818 |
< |
lennard-Jones sphere solvents. The solvent molecules are excluded from |
| 819 |
< |
the system, the experimental value of water viscosity is applied to |
| 820 |
< |
mimic the heat bath. Fig. XXX is the snapshot of the stable |
| 821 |
< |
configuration of the system, the ripple structure stayed stable after |
| 822 |
< |
100 ns run. The efficiency of the simulation is increased by one order |
| 823 |
< |
of magnitude. |
| 1273 |
> |
\section{Application: A rigid-body lipid bilayer} |
| 1274 |
|
|
| 1275 |
< |
\subsection{Langevin Dynamics of Banana Shaped Molecules} |
| 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 |
|
|
| 827 |
– |
In order to verify that Langevin dynamics can mimic the dynamics of |
| 828 |
– |
the systems absent of explicit solvents, we carried out two sets of |
| 829 |
– |
simulations and compare their dynamic properties. |
| 830 |
– |
Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation |
| 831 |
– |
made of 256 pentane molecules and two banana shaped molecules at |
| 832 |
– |
273~K. It has an equivalent implicit solvent system containing only |
| 833 |
– |
two banana shaped molecules with viscosity of 0.289 center poise. To |
| 834 |
– |
calculate the hydrodynamic properties of the banana shaped molecule, |
| 835 |
– |
we created a rough shell model (see Fig.~\ref{langevin:roughShell}), |
| 836 |
– |
in which the banana shaped molecule is represented as a ``shell'' |
| 837 |
– |
made of 2266 small identical beads with size of 0.3 \AA on the |
| 838 |
– |
surface. Applying the procedure described in |
| 839 |
– |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
| 840 |
– |
identified the center of resistance at (0 $\rm{\AA}$, 0.7482 $\rm{\AA}$, |
| 841 |
– |
-0.1988 $\rm{\AA}$), as well as the resistance tensor, |
| 842 |
– |
\[ |
| 843 |
– |
\left( {\begin{array}{*{20}c} |
| 844 |
– |
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
| 845 |
– |
0& 0.9270&-0.007063& 0.08585&0&0\\ |
| 846 |
– |
0&-0.007063&0.7494&0.2057&0&0\\ |
| 847 |
– |
0&0.0858&0.2057& 58.64& 0&0\\ |
| 848 |
– |
0.08585&0&0&0&48.30&3.219&\\ |
| 849 |
– |
0.2057&0&0&0&3.219&10.7373\\ |
| 850 |
– |
\end{array}} \right). |
| 851 |
– |
\] |
| 852 |
– |
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. |
| 853 |
– |
Curves of the velocity auto-correlation functions in |
| 854 |
– |
Fig.~\ref{langevin:vacf} were shown to match each other very well. |
| 855 |
– |
However, because of the stochastic nature, simulation using Langevin |
| 856 |
– |
dynamics was shown to decay slightly faster than MD. In order to |
| 857 |
– |
study the rotational motion of the molecules, we also calculated the |
| 858 |
– |
auto-correlation function of the principle axis of the second GB |
| 859 |
– |
particle, $u$. The discrepancy shown in Fig.~\ref{langevin:uacf} was |
| 860 |
– |
probably due to the reason that we used the experimental viscosity directly instead of calculating bulk viscosity from simulation. |
| 861 |
– |
|
| 1299 |
|
\begin{figure} |
| 1300 |
|
\centering |
| 1301 |
< |
\includegraphics[width=\linewidth]{roughShell.pdf} |
| 1302 |
< |
\caption[Rough shell model for banana shaped molecule]{Rough shell |
| 1303 |
< |
model for banana shaped molecule.} \label{langevin:roughShell} |
| 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 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 |
|
|
| 869 |
– |
\begin{figure} |
| 870 |
– |
\centering |
| 871 |
– |
\includegraphics[width=\linewidth]{twoBanana.pdf} |
| 872 |
– |
\caption[Snapshot from Simulation of Two Banana Shaped Molecules and |
| 873 |
– |
256 Pentane Molecules]{Snapshot from simulation of two Banana shaped |
| 874 |
– |
molecules and 256 pentane molecules.} \label{langevin:twoBanana} |
| 875 |
– |
\end{figure} |
| 876 |
– |
|
| 877 |
– |
\begin{figure} |
| 878 |
– |
\centering |
| 879 |
– |
\includegraphics[width=\linewidth]{vacf.pdf} |
| 880 |
– |
\caption[Plots of Velocity Auto-correlation Functions]{Velocity |
| 881 |
– |
auto-correlation functions of NVE (explicit solvent) in blue and |
| 882 |
– |
Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf} |
| 883 |
– |
\end{figure} |
| 884 |
– |
|
| 885 |
– |
\begin{figure} |
| 886 |
– |
\centering |
| 887 |
– |
\includegraphics[width=\linewidth]{uacf.pdf} |
| 888 |
– |
\caption[Auto-correlation functions of the principle axis of the |
| 889 |
– |
middle GB particle]{Auto-correlation functions of the principle axis |
| 890 |
– |
of the middle GB particle of NVE (blue) and Langevin dynamics |
| 891 |
– |
(red).} \label{langevin:uacf} |
| 892 |
– |
\end{figure} |
| 893 |
– |
|
| 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. The temperature control |
| 1313 |
< |
ability of this algorithm was demonstrated by a set of simulations |
| 1314 |
< |
with different viscosities. It was also shown to have significant |
| 1315 |
< |
advantage of producing rapid thermal equilibration over |
| 1316 |
< |
Nos\'{e}-Hoover method. Further studies in systems involving banana |
| 1317 |
< |
shaped molecules illustrated that the dynamic properties could be |
| 1318 |
< |
preserved by using this new algorithm as an 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{jcp2} |
| 1356 |
|
\bibliography{langevin} |
| 916 |
– |
|
| 1357 |
|
\end{document} |
