4 |
|
\usepackage{endfloat} |
5 |
|
\usepackage{amsmath,bm} |
6 |
|
\usepackage{amssymb} |
7 |
– |
\usepackage{epsf} |
7 |
|
\usepackage{times} |
8 |
|
\usepackage{mathptmx} |
9 |
|
\usepackage{setspace} |
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 |
20 |
> |
\renewcommand\citemid{\ } % no comma in optional referenc note |
21 |
|
|
22 |
|
\begin{document} |
23 |
|
|
24 |
< |
\title{Langevin Dynamics for Rigid Body of Arbitrary Shape } |
24 |
> |
\title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape } |
25 |
|
|
26 |
< |
\author{Teng Lin, Xiuquan Sun and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
26 |
> |
\author{Xiuquan Sun, Teng Lin and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
27 |
|
gezelter@nd.edu} \\ |
28 |
|
Department of Chemistry and Biochemistry\\ |
29 |
|
University of Notre Dame\\ |
170 |
|
algorithm for arbitrary-shaped rigid particles by integrating the |
171 |
|
accurate estimation of friction tensor from hydrodynamics theory |
172 |
|
into the sophisticated rigid body dynamics algorithms. |
174 |
– |
|
175 |
– |
\section{Computational Methods{\label{methodSec}}} |
173 |
|
|
174 |
|
\subsection{\label{introSection:frictionTensor}Friction Tensor} |
175 |
|
Theoretically, the friction kernel can be determined using the |
365 |
|
\begin{eqnarray} |
366 |
|
\Xi _{}^{tt} & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\ |
367 |
|
\Xi _{}^{tr} & = & \Xi _{}^{rt} = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\ |
368 |
< |
\Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\ |
368 |
> |
\Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } |
369 |
> |
U_j + 6 \eta V {\bf I}. \notag |
370 |
|
\label{introEquation:ResistanceTensorArbitraryOrigin} |
371 |
|
\end{eqnarray} |
372 |
+ |
The final term in the expression for $\Xi^{rr}$ is correction that |
373 |
+ |
accounts for errors in the rotational motion of certain kinds of bead |
374 |
+ |
models. The additive correction uses the solvent viscosity ($\eta$) |
375 |
+ |
as well as the total volume of the beads that contribute to the |
376 |
+ |
hydrodynamic model, |
377 |
+ |
\begin{equation} |
378 |
+ |
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3, |
379 |
+ |
\end{equation} |
380 |
+ |
where $\sigma_i$ is the radius of bead $i$. This correction term was |
381 |
+ |
rigorously tested and compared with the analytical results for |
382 |
+ |
two-sphere and ellipsoidal systems by Garcia de la Torre and |
383 |
+ |
Rodes.\cite{Torre:1983lr} |
384 |
+ |
|
385 |
+ |
|
386 |
|
The resistance tensor depends on the origin to which they refer. The |
387 |
|
proper location for applying the friction force is the center of |
388 |
|
resistance (or center of reaction), at which the trace of rotational |
439 |
|
where $x_OR$, $y_OR$, $z_OR$ are the components of the vector |
440 |
|
joining center of resistance $R$ and origin $O$. |
441 |
|
|
430 |
– |
\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
442 |
|
|
443 |
+ |
\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} |
444 |
|
Consider the Langevin equations of motion in generalized coordinates |
445 |
|
\begin{equation} |
446 |
|
M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)} + F_{r,i} (t) |
587 |
|
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
588 |
|
\end{align*} |
589 |
|
|
590 |
< |
\section{Results and Discussion} |
591 |
< |
|
592 |
< |
The Langevin algorithm described in previous section has been |
593 |
< |
implemented in {\sc oopse}\cite{Meineke2005} and applied to studies |
594 |
< |
of the static and dynamic properties in several systems. |
595 |
< |
|
596 |
< |
\subsection{Temperature Control} |
590 |
> |
\section{Validating the Method\label{sec:validating}} |
591 |
> |
In order to validate our Langevin integrator for arbitrarily-shaped |
592 |
> |
rigid bodies, we implemented the algorithm in {\sc |
593 |
> |
oopse}\cite{Meineke2005} and compared the results of this algorithm |
594 |
> |
with the known |
595 |
> |
hydrodynamic limiting behavior for a few model systems, and to |
596 |
> |
microcanonical molecular dynamics simulations for some more |
597 |
> |
complicated bodies. The model systems and their analytical behavior |
598 |
> |
(if known) are summarized below. Parameters for the primary particles |
599 |
> |
comprising our model systems are given in table \ref{tab:parameters}, |
600 |
> |
and a sketch of the arrangement of these primary particles into the |
601 |
> |
model rigid bodies is shown in figure \ref{fig:models}. In table |
602 |
> |
\ref{tab:parameters}, $d$ and $l$ are the physical dimensions of |
603 |
> |
ellipsoidal (Gay-Berne) particles. For spherical particles, the value |
604 |
> |
of the Lennard-Jones $\sigma$ parameter is the particle diameter |
605 |
> |
($d$). Gay-Berne ellipsoids have an energy scaling parameter, |
606 |
> |
$\epsilon^s$, which describes the well depth for two identical |
607 |
> |
ellipsoids in a {\it side-by-side} configuration. Additionally, a |
608 |
> |
well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, |
609 |
> |
describes the ratio between the well depths in the {\it end-to-end} |
610 |
> |
and side-by-side configurations. For spheres, $\epsilon^r \equiv 1$. |
611 |
> |
Moments of inertia are also required to describe the motion of primary |
612 |
> |
particles with orientational degrees of freedom. |
613 |
|
|
614 |
< |
As shown in Eq.~\ref{randomForce}, random collisions associated with |
615 |
< |
the solvent's thermal motions is controlled by the external |
588 |
< |
temperature. The capability to maintain the temperature of the whole |
589 |
< |
system was usually used to measure the stability and efficiency of |
590 |
< |
the algorithm. In order to verify the stability of this new |
591 |
< |
algorithm, a series of simulations are performed on system |
592 |
< |
consisiting of 256 SSD water molecules with different viscosities. |
593 |
< |
The initial configuration for the simulations is taken from a 1ns |
594 |
< |
NVT simulation with a cubic box of 19.7166~\AA. All simulation are |
595 |
< |
carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns |
596 |
< |
with reference temperature at 300~K. The average temperature as a |
597 |
< |
function of $\eta$ is shown in Table \ref{langevin:viscosity} where |
598 |
< |
the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 - |
599 |
< |
1$ poise. The better temperature control at higher viscosity can be |
600 |
< |
explained by the finite size effect and relative slow relaxation |
601 |
< |
rate at lower viscosity regime. |
602 |
< |
\begin{table} |
603 |
< |
\caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF |
604 |
< |
SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.} |
605 |
< |
\label{langevin:viscosity} |
614 |
> |
\begin{table*} |
615 |
> |
\begin{minipage}{\linewidth} |
616 |
|
\begin{center} |
617 |
< |
\begin{tabular}{lll} |
618 |
< |
\hline |
619 |
< |
$\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\ |
620 |
< |
\hline |
621 |
< |
1 & 300.47 & 10.99 \\ |
622 |
< |
0.1 & 301.19 & 11.136 \\ |
623 |
< |
0.01 & 303.04 & 11.796 \\ |
624 |
< |
\hline |
617 |
> |
\caption{Parameters for the primary particles in use by the rigid body |
618 |
> |
models in figure \ref{fig:models}.} |
619 |
> |
\begin{tabular}{lrcccccccc} |
620 |
> |
\hline |
621 |
> |
& & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\ |
622 |
> |
& & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ (kcal/mol) & $\epsilon^r$ & |
623 |
> |
$m$ (amu) & $I_{xx}$ & $I_{yy}$ & $I_{zz}$ \\ \hline |
624 |
> |
Sphere & & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
625 |
> |
Ellipsoid & & 4.6 & 13.8 & 0.8 & 0.2 & 200 & 2105 & 2105 & 421 \\ |
626 |
> |
Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
627 |
> |
Banana &(3 identical ellipsoids)& 4.2 & 11.2 & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\ |
628 |
> |
Lipid: & Spherical Head & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\ |
629 |
> |
& Ellipsoidal Tail & 4.6 & 13.8 & 0.8 & 0.2 & 760 & 45000 & 45000 & 9000 \\ |
630 |
> |
Solvent & & 4.7 & $= d$ & 0.8 & 1 & 72.06 & & & \\ |
631 |
> |
\hline |
632 |
|
\end{tabular} |
633 |
+ |
\label{tab:parameters} |
634 |
|
\end{center} |
635 |
< |
\end{table} |
635 |
> |
\end{minipage} |
636 |
> |
\end{table*} |
637 |
|
|
619 |
– |
Another set of calculations were performed to study the efficiency of |
620 |
– |
temperature control using different temperature coupling schemes. |
621 |
– |
The starting configuration is cooled to 173~K and evolved using NVE, |
622 |
– |
NVT, and Langevin dynamic with time step of 2 fs. |
623 |
– |
Fig.~\ref{langevin:temperature} shows the heating curve obtained as |
624 |
– |
the systems reach equilibrium. The orange curve in |
625 |
– |
Fig.~\ref{langevin:temperature} represents the simulation using |
626 |
– |
Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps |
627 |
– |
which gives reasonable tight coupling, while the blue one from |
628 |
– |
Langevin dynamics with viscosity of 0.1 poise demonstrates a faster |
629 |
– |
scaling to the desire temperature. When $ \eta = 0$, Langevin dynamics becomes normal |
630 |
– |
NVE (see orange curve in Fig.~\ref{langevin:temperature}) which |
631 |
– |
loses the temperature control ability. |
632 |
– |
|
638 |
|
\begin{figure} |
639 |
|
\centering |
640 |
< |
\includegraphics[width=\linewidth]{temperature.eps} |
641 |
< |
\caption[Plot of Temperature Fluctuation Versus Time]{Plot of |
642 |
< |
temperature fluctuation versus time.} \label{langevin:temperature} |
640 |
> |
\includegraphics[width=3in]{sketch} |
641 |
> |
\caption[Sketch of the model systems]{A sketch of the model systems |
642 |
> |
used in evaluating the behavior of the rigid body Langevin |
643 |
> |
integrator.} \label{fig:models} |
644 |
|
\end{figure} |
645 |
|
|
646 |
< |
\subsection{Langevin Dynamics of Banana Shaped Molecules} |
646 |
> |
\subsection{Simulation Methodology} |
647 |
> |
We performed reference microcanonical simulations with explicit |
648 |
> |
solvents for each of the different model system. In each case there |
649 |
> |
was one solute model and 1929 solvent molecules present in the |
650 |
> |
simulation box. All simulations were equilibrated using a |
651 |
> |
constant-pressure and temperature integrator with target values of 300 |
652 |
> |
K for the temperature and 1 atm for pressure. Following this stage, |
653 |
> |
further equilibration and sampling was done in a microcanonical |
654 |
> |
ensemble. Since the model bodies are typically quite massive, we were |
655 |
> |
able to use a time step of 25 fs. |
656 |
|
|
657 |
< |
In order to verify that Langevin dynamics can mimic the dynamics of |
658 |
< |
the systems absent of explicit solvents, we carried out two sets of |
659 |
< |
simulations and compare their dynamic properties. |
660 |
< |
Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation |
661 |
< |
made of 256 pentane molecules and two banana shaped molecules at |
662 |
< |
273~K. It has an equivalent implicit solvent system containing only |
663 |
< |
two banana shaped molecules with viscosity of 0.289 center poise. To |
664 |
< |
calculate the hydrodynamic properties of the banana shaped molecule, |
665 |
< |
we created a rough shell model (see Fig.~\ref{langevin:roughShell}), |
666 |
< |
in which the banana shaped molecule is represented as a ``shell'' |
667 |
< |
made of 2266 small identical beads with size of 0.3 \AA on the |
668 |
< |
surface. Applying the procedure described in |
669 |
< |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
670 |
< |
identified the center of resistance at (0 $\rm{\AA}$, 0.7482 $\rm{\AA}$, |
671 |
< |
-0.1988 $\rm{\AA}$), as well as the resistance tensor, |
672 |
< |
\[ |
673 |
< |
\left( {\begin{array}{*{20}c} |
674 |
< |
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
660 |
< |
0& 0.9270&-0.007063& 0.08585&0&0\\ |
661 |
< |
0&-0.007063&0.7494&0.2057&0&0\\ |
662 |
< |
0&0.0858&0.2057& 58.64& 0&0\\ |
663 |
< |
0.08585&0&0&0&48.30&3.219&\\ |
664 |
< |
0.2057&0&0&0&3.219&10.7373\\ |
665 |
< |
\end{array}} \right). |
666 |
< |
\] |
667 |
< |
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. |
668 |
< |
Curves of the velocity auto-correlation functions in |
669 |
< |
Fig.~\ref{langevin:vacf} were shown to match each other very well. |
670 |
< |
However, because of the stochastic nature, simulation using Langevin |
671 |
< |
dynamics was shown to decay slightly faster than MD. In order to |
672 |
< |
study the rotational motion of the molecules, we also calculated the |
673 |
< |
auto-correlation function of the principle axis of the second GB |
674 |
< |
particle, $u$. The discrepancy shown in Fig.~\ref{langevin:uacf} was |
675 |
< |
probably due to the reason that we used the experimental viscosity directly instead of calculating bulk viscosity from simulation. |
657 |
> |
The model systems studied used both Lennard-Jones spheres as well as |
658 |
> |
uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne |
659 |
> |
potential was a single site model for the interactions of rigid |
660 |
> |
ellipsoidal molecules.\cite{Gay81} It can be thought of as a |
661 |
> |
modification of the Gaussian overlap model originally described by |
662 |
> |
Berne and Pechukas.\cite{Berne72} The potential is constructed in the |
663 |
> |
familiar form of the Lennard-Jones function using |
664 |
> |
orientation-dependent $\sigma$ and $\epsilon$ parameters, |
665 |
> |
\begin{equation*} |
666 |
> |
V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat |
667 |
> |
r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, |
668 |
> |
{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u |
669 |
> |
}_i}, |
670 |
> |
{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12} |
671 |
> |
-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, |
672 |
> |
{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right] |
673 |
> |
\label{eq:gb} |
674 |
> |
\end{equation*} |
675 |
|
|
676 |
< |
\begin{figure} |
677 |
< |
\centering |
678 |
< |
\includegraphics[width=\linewidth]{roughShell.eps} |
679 |
< |
\caption[Rough shell model for banana shaped molecule]{Rough shell |
680 |
< |
model for banana shaped molecule.} \label{langevin:roughShell} |
681 |
< |
\end{figure} |
676 |
> |
The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf |
677 |
> |
\hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf |
678 |
> |
\hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters |
679 |
> |
are dependent on the relative orientations of the two ellipsoids (${\bf |
680 |
> |
\hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the |
681 |
> |
inter-ellipsoid separation (${\bf \hat{r}}_{ij}$). The shape and |
682 |
> |
attractiveness of each ellipsoid is governed by a relatively small set |
683 |
> |
of parameters: $l$ and $d$ describe the length and width of each |
684 |
> |
uniaxial ellipsoid, while $\epsilon^s$, which describes the well depth |
685 |
> |
for two identical ellipsoids in a {\it side-by-side} configuration. |
686 |
> |
Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e / |
687 |
> |
\epsilon^s$, describes the ratio between the well depths in the {\it |
688 |
> |
end-to-end} and side-by-side configurations. Details of the potential |
689 |
> |
are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an |
690 |
> |
excellent overview of the computational methods that can be used to |
691 |
> |
efficiently compute forces and torques for this potential can be found |
692 |
> |
in Ref. \citen{Golubkov06} |
693 |
> |
|
694 |
> |
For the interaction between nonequivalent uniaxial ellipsoids (or |
695 |
> |
between spheres and ellipsoids), the spheres are treated as ellipsoids |
696 |
> |
with an aspect ratio of 1 ($d = l$) and with an well depth ratio |
697 |
> |
($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$). The form of the |
698 |
> |
Gay-Berne potential we are using was generalized by Cleaver {\it et |
699 |
> |
al.} and is appropriate for dissimilar uniaxial |
700 |
> |
ellipsoids.\cite{Cleaver96} |
701 |
> |
|
702 |
> |
A switching function was applied to all potentials to smoothly turn |
703 |
> |
off the interactions between a range of $22$ and $25$ \AA. The |
704 |
> |
switching function was the standard (cubic) function, |
705 |
> |
\begin{equation} |
706 |
> |
s(r) = |
707 |
> |
\begin{cases} |
708 |
> |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
709 |
> |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
710 |
> |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
711 |
> |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
712 |
> |
0 & \text{if $r > r_{\text{cut}}$.} |
713 |
> |
\end{cases} |
714 |
> |
\label{eq:switchingFunc} |
715 |
> |
\end{equation} |
716 |
> |
|
717 |
> |
To measure shear viscosities from our microcanonical simulations, we |
718 |
> |
used the Einstein form of the pressure correlation function,\cite{hess:209} |
719 |
> |
\begin{equation} |
720 |
> |
\eta = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left( |
721 |
> |
\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}. |
722 |
> |
\label{eq:shear} |
723 |
> |
\end{equation} |
724 |
> |
A similar form exists for the bulk viscosity |
725 |
> |
\begin{equation} |
726 |
> |
\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left( |
727 |
> |
\int_{t_0}^{t_0 + t} |
728 |
> |
\left(P\left(t'\right)-\left\langle P \right\rangle \right)dt' |
729 |
> |
\right)^2 \right\rangle_{t_0}. |
730 |
> |
\end{equation} |
731 |
> |
Alternatively, the shear viscosity can also be calculated using a |
732 |
> |
Green-Kubo formula with the off-diagonal pressure tensor correlation function, |
733 |
> |
\begin{equation} |
734 |
> |
\eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0 |
735 |
> |
+ t) \right\rangle_{t_0} dt, |
736 |
> |
\end{equation} |
737 |
> |
although this method converges extremely slowly and is not practical |
738 |
> |
for obtaining viscosities from molecular dynamics simulations. |
739 |
> |
|
740 |
> |
The Langevin dynamics for the different model systems were performed |
741 |
> |
at the same temperature as the average temperature of the |
742 |
> |
microcanonical simulations and with a solvent viscosity taken from |
743 |
> |
Eq. (\ref{eq:shear}) applied to these simulations. We used 1024 |
744 |
> |
independent solute simulations to obtain statistics on our Langevin |
745 |
> |
integrator. |
746 |
> |
|
747 |
> |
\subsection{Analysis} |
748 |
> |
|
749 |
> |
The quantities of interest when comparing the Langevin integrator to |
750 |
> |
analytic hydrodynamic equations and to molecular dynamics simulations |
751 |
> |
are typically translational diffusion constants and orientational |
752 |
> |
relaxation times. Translational diffusion constants for point |
753 |
> |
particles are computed easily from the long-time slope of the |
754 |
> |
mean-square displacement, |
755 |
> |
\begin{equation} |
756 |
> |
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \left\langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \right\rangle, |
757 |
> |
\end{equation} |
758 |
> |
of the solute molecules. For models in which the translational |
759 |
> |
diffusion tensor (${\bf D}_{tt}$) has non-degenerate eigenvalues |
760 |
> |
(i.e. any non-spherically-symmetric rigid body), it is possible to |
761 |
> |
compute the diffusive behavior for motion parallel to each body-fixed |
762 |
> |
axis by projecting the displacement of the particle onto the |
763 |
> |
body-fixed reference frame at $t=0$. With an isotropic solvent, as we |
764 |
> |
have used in this study, there are differences between the three |
765 |
> |
diffusion constants, but these must converge to the same value at |
766 |
> |
longer times. Translational diffusion constants for the different |
767 |
> |
shaped models are shown in table \ref{tab:translation}. |
768 |
> |
|
769 |
> |
In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational |
770 |
> |
diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object |
771 |
> |
{\it around} a particular body-fixed axis and {\it not} the diffusion |
772 |
> |
of a vector pointing along the axis. However, these eigenvalues can |
773 |
> |
be combined to find 5 characteristic rotational relaxation |
774 |
> |
times,\cite{PhysRev.119.53,Berne90} |
775 |
> |
\begin{eqnarray} |
776 |
> |
1 / \tau_1 & = & 6 D_r + 2 \Delta \\ |
777 |
> |
1 / \tau_2 & = & 6 D_r - 2 \Delta \\ |
778 |
> |
1 / \tau_3 & = & 3 (D_r + D_1) \\ |
779 |
> |
1 / \tau_4 & = & 3 (D_r + D_2) \\ |
780 |
> |
1 / \tau_5 & = & 3 (D_r + D_3) |
781 |
> |
\end{eqnarray} |
782 |
> |
where |
783 |
> |
\begin{equation} |
784 |
> |
D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right) |
785 |
> |
\end{equation} |
786 |
> |
and |
787 |
> |
\begin{equation} |
788 |
> |
\Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2} |
789 |
> |
\end{equation} |
790 |
> |
Each of these characteristic times can be used to predict the decay of |
791 |
> |
part of the rotational correlation function when $\ell = 2$, |
792 |
> |
\begin{equation} |
793 |
> |
C_2(t) = \frac{a^2}{N^2} e^{-t/\tau_1} + \frac{b^2}{N^2} e^{-t/\tau_2}. |
794 |
> |
\end{equation} |
795 |
> |
This is the same as the $F^2_{0,0}(t)$ correlation function that |
796 |
> |
appears in Ref. \citen{Berne90}. The amplitudes of the two decay |
797 |
> |
terms are expressed in terms of three dimensionless functions of the |
798 |
> |
eigenvalues: $a = \sqrt{3} (D_1 - D_2)$, $b = (2D_3 - D_1 - D_2 + |
799 |
> |
2\Delta)$, and $N = 2 \sqrt{\Delta b}$. Similar expressions can be |
800 |
> |
obtained for other angular momentum correlation |
801 |
> |
functions.\cite{PhysRev.119.53,Berne90} In all of the model systems we |
802 |
> |
studied, only one of the amplitudes of the two decay terms was |
803 |
> |
non-zero, so it was possible to derive a single relaxation time for |
804 |
> |
each of the hydrodynamic tensors. In many cases, these characteristic |
805 |
> |
times are averaged and reported in the literature as a single relaxation |
806 |
> |
time,\cite{Garcia-de-la-Torre:1997qy} |
807 |
> |
\begin{equation} |
808 |
> |
1 / \tau_0 = \frac{1}{5} \sum_{i=1}^5 \tau_{i}^{-1}, |
809 |
> |
\end{equation} |
810 |
> |
although for the cases reported here, this averaging is not necessary |
811 |
> |
and only one of the five relaxation times is relevant. |
812 |
> |
|
813 |
> |
To test the Langevin integrator's behavior for rotational relaxation, |
814 |
> |
we have compared the analytical orientational relaxation times (if |
815 |
> |
they are known) with the general result from the diffusion tensor and |
816 |
> |
with the results from both the explicitly solvated molecular dynamics |
817 |
> |
and Langevin simulations. Relaxation times from simulations (both |
818 |
> |
microcanonical and Langevin), were computed using Legendre polynomial |
819 |
> |
correlation functions for a unit vector (${\bf u}$) fixed along one or |
820 |
> |
more of the body-fixed axes of the model. |
821 |
> |
\begin{equation} |
822 |
> |
C_{\ell}(t) = \left\langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf |
823 |
> |
u}_{i}(0) \right) \right\rangle |
824 |
> |
\end{equation} |
825 |
> |
For simulations in the high-friction limit, orientational correlation |
826 |
> |
times can then be obtained from exponential fits of this function, or by |
827 |
> |
integrating, |
828 |
> |
\begin{equation} |
829 |
> |
\tau = \ell (\ell + 1) \int_0^{\infty} C_{\ell}(t) dt. |
830 |
> |
\end{equation} |
831 |
> |
In lower-friction solvents, the Legendre correlation functions often |
832 |
> |
exhibit non-exponential decay, and may not be characterized by a |
833 |
> |
single decay constant. |
834 |
> |
|
835 |
> |
In table \ref{tab:rotation} we show the characteristic rotational |
836 |
> |
relaxation times (based on the diffusion tensor) for each of the model |
837 |
> |
systems compared with the values obtained via microcanonical and Langevin |
838 |
> |
simulations. |
839 |
|
|
840 |
+ |
\subsection{Spherical particles} |
841 |
+ |
Our model system for spherical particles was a Lennard-Jones sphere of |
842 |
+ |
diameter ($\sigma$) 6.5 \AA\ in a sea of smaller spheres ($\sigma$ = |
843 |
+ |
4.7 \AA). The well depth ($\epsilon$) for both particles was set to |
844 |
+ |
an arbitrary value of 0.8 kcal/mol. |
845 |
+ |
|
846 |
+ |
The Stokes-Einstein behavior of large spherical particles in |
847 |
+ |
hydrodynamic flows is well known, giving translational friction |
848 |
+ |
coefficients of $6 \pi \eta R$ (stick boundary conditions) and |
849 |
+ |
rotational friction coefficients of $8 \pi \eta R^3$. Recently, |
850 |
+ |
Schmidt and Skinner have computed the behavior of spherical tag |
851 |
+ |
particles in molecular dynamics simulations, and have shown that {\it |
852 |
+ |
slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more |
853 |
+ |
appropriate for molecule-sized spheres embedded in a sea of spherical |
854 |
+ |
solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
855 |
+ |
|
856 |
+ |
Our simulation results show similar behavior to the behavior observed |
857 |
+ |
by Schmidt and Skinner. The diffusion constant obtained from our |
858 |
+ |
microcanonical molecular dynamics simulations lies between the slip |
859 |
+ |
and stick boundary condition results obtained via Stokes-Einstein |
860 |
+ |
behavior. Since the Langevin integrator assumes Stokes-Einstein stick |
861 |
+ |
boundary conditions in calculating the drag and random forces for |
862 |
+ |
spherical particles, our Langevin routine obtains nearly quantitative |
863 |
+ |
agreement with the hydrodynamic results for spherical particles. One |
864 |
+ |
avenue for improvement of the method would be to compute elements of |
865 |
+ |
$\Xi_{tt}$ assuming behavior intermediate between the two boundary |
866 |
+ |
conditions. |
867 |
+ |
|
868 |
+ |
In the explicit solvent simulations, both our solute and solvent |
869 |
+ |
particles were structureless, exerting no torques upon each other. |
870 |
+ |
Therefore, there are not rotational correlation times available for |
871 |
+ |
this model system. |
872 |
+ |
|
873 |
+ |
\subsection{Ellipsoids} |
874 |
+ |
For uniaxial ellipsoids ($a > b = c$), Perrin's formulae for both |
875 |
+ |
translational and rotational diffusion of each of the body-fixed axes |
876 |
+ |
can be combined to give a single translational diffusion |
877 |
+ |
constant,\cite{Berne90} |
878 |
+ |
\begin{equation} |
879 |
+ |
D = \frac{k_B T}{6 \pi \eta a} G(\rho), |
880 |
+ |
\label{Dperrin} |
881 |
+ |
\end{equation} |
882 |
+ |
as well as a single rotational diffusion coefficient, |
883 |
+ |
\begin{equation} |
884 |
+ |
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2) |
885 |
+ |
G(\rho) - 1}{1 - \rho^4} \right\}. |
886 |
+ |
\label{ThetaPerrin} |
887 |
+ |
\end{equation} |
888 |
+ |
In these expressions, $G(\rho)$ is a function of the axial ratio |
889 |
+ |
($\rho = b / a$), which for prolate ellipsoids, is |
890 |
+ |
\begin{equation} |
891 |
+ |
G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 - |
892 |
+ |
\rho^2)^{1/2}}{\rho} \right\} |
893 |
+ |
\label{GPerrin} |
894 |
+ |
\end{equation} |
895 |
+ |
Again, there is some uncertainty about the correct boundary conditions |
896 |
+ |
to use for molecular-scale ellipsoids in a sea of similarly-sized |
897 |
+ |
solvent particles. Ravichandran and Bagchi found that {\it slip} |
898 |
+ |
boundary conditions most closely resembled the simulation |
899 |
+ |
results,\cite{Ravichandran:1999fk} in agreement with earlier work of |
900 |
+ |
Tang and Evans.\cite{TANG:1993lr} |
901 |
+ |
|
902 |
+ |
Even though there are analytic resistance tensors for ellipsoids, we |
903 |
+ |
constructed a rough-shell model using 2135 beads (each with a diameter |
904 |
+ |
of 0.25 \AA) to approximate the shape of the model ellipsoid. We |
905 |
+ |
compared the Langevin dynamics from both the simple ellipsoidal |
906 |
+ |
resistance tensor and the rough shell approximation with |
907 |
+ |
microcanonical simulations and the predictions of Perrin. As in the |
908 |
+ |
case of our spherical model system, the Langevin integrator reproduces |
909 |
+ |
almost exactly the behavior of the Perrin formulae (which is |
910 |
+ |
unsurprising given that the Perrin formulae were used to derive the |
911 |
+ |
drag and random forces applied to the ellipsoid). We obtain |
912 |
+ |
translational diffusion constants and rotational correlation times |
913 |
+ |
that are within a few percent of the analytic values for both the |
914 |
+ |
exact treatment of the diffusion tensor as well as the rough-shell |
915 |
+ |
model for the ellipsoid. |
916 |
+ |
|
917 |
+ |
The translational diffusion constants from the microcanonical simulations |
918 |
+ |
agree well with the predictions of the Perrin model, although the rotational |
919 |
+ |
correlation times are a factor of 2 shorter than expected from hydrodynamic |
920 |
+ |
theory. One explanation for the slower rotation |
921 |
+ |
of explicitly-solvated ellipsoids is the possibility that solute-solvent |
922 |
+ |
collisions happen at both ends of the solute whenever the principal |
923 |
+ |
axis of the ellipsoid is turning. In the upper portion of figure |
924 |
+ |
\ref{fig:explanation} we sketch a physical picture of this explanation. |
925 |
+ |
Since our Langevin integrator is providing nearly quantitative agreement with |
926 |
+ |
the Perrin model, it also predicts orientational diffusion for ellipsoids that |
927 |
+ |
exceed explicitly solvated correlation times by a factor of two. |
928 |
+ |
|
929 |
+ |
\subsection{Rigid dumbbells} |
930 |
+ |
Perhaps the only {\it composite} rigid body for which analytic |
931 |
+ |
expressions for the hydrodynamic tensor are available is the |
932 |
+ |
two-sphere dumbbell model. This model consists of two non-overlapping |
933 |
+ |
spheres held by a rigid bond connecting their centers. There are |
934 |
+ |
competing expressions for the 6x6 resistance tensor for this |
935 |
+ |
model. Equation (\ref{introEquation:oseenTensor}) above gives the |
936 |
+ |
original Oseen tensor, while the second order expression introduced by |
937 |
+ |
Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la |
938 |
+ |
Torre and Bloomfield,\cite{Torre1977} is given above as |
939 |
+ |
Eq. (\ref{introEquation:RPTensorNonOverlapped}). In our case, we use |
940 |
+ |
a model dumbbell in which the two spheres are identical Lennard-Jones |
941 |
+ |
particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at |
942 |
+ |
a distance of 6.532 \AA. |
943 |
+ |
|
944 |
+ |
The theoretical values for the translational diffusion constant of the |
945 |
+ |
dumbbell are calculated from the work of Stimson and Jeffery, who |
946 |
+ |
studied the motion of this system in a flow parallel to the |
947 |
+ |
inter-sphere axis,\cite{Stimson:1926qy} and Davis, who studied the |
948 |
+ |
motion in a flow {\it perpendicular} to the inter-sphere |
949 |
+ |
axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it |
950 |
+ |
orientational} correlation times for this model system (other than |
951 |
+ |
those derived from the 6 x 6 tensors mentioned above). |
952 |
+ |
|
953 |
+ |
The bead model for this model system comprises the two large spheres |
954 |
+ |
by themselves, while the rough shell approximation used 3368 separate |
955 |
+ |
beads (each with a diameter of 0.25 \AA) to approximate the shape of |
956 |
+ |
the rigid body. The hydrodynamics tensors computed from both the bead |
957 |
+ |
and rough shell models are remarkably similar. Computing the initial |
958 |
+ |
hydrodynamic tensor for a rough shell model can be quite expensive (in |
959 |
+ |
this case it requires inverting a 10104 x 10104 matrix), while the |
960 |
+ |
bead model is typically easy to compute (in this case requiring |
961 |
+ |
inversion of a 6 x 6 matrix). |
962 |
+ |
|
963 |
|
\begin{figure} |
964 |
|
\centering |
965 |
< |
\includegraphics[width=\linewidth]{twoBanana.eps} |
966 |
< |
\caption[Snapshot from Simulation of Two Banana Shaped Molecules and |
967 |
< |
256 Pentane Molecules]{Snapshot from simulation of two Banana shaped |
968 |
< |
molecules and 256 pentane molecules.} \label{langevin:twoBanana} |
965 |
> |
\includegraphics[width=2in]{RoughShell} |
966 |
> |
\caption[Model rigid bodies and their rough shell approximations]{The |
967 |
> |
model rigid bodies (left column) used to test this algorithm and their |
968 |
> |
rough-shell approximations (right-column) that were used to compute |
969 |
> |
the hydrodynamic tensors. The top two models (ellipsoid and dumbbell) |
970 |
> |
have analytic solutions and were used to test the rough shell |
971 |
> |
approximation. The lower two models (banana and lipid) were compared |
972 |
> |
with explicitly-solvated molecular dynamics simulations. } |
973 |
> |
\label{fig:roughShell} |
974 |
|
\end{figure} |
975 |
|
|
976 |
+ |
|
977 |
+ |
Once the hydrodynamic tensor has been computed, there is no additional |
978 |
+ |
penalty for carrying out a Langevin simulation with either of the two |
979 |
+ |
different hydrodynamics models. Our naive expectation is that since |
980 |
+ |
the rigid body's surface is roughened under the various shell models, |
981 |
+ |
the diffusion constants will be even farther from the ``slip'' |
982 |
+ |
boundary conditions than observed for the bead model (which uses a |
983 |
+ |
Stokes-Einstein model to arrive at the hydrodynamic tensor). For the |
984 |
+ |
dumbbell, this prediction is correct although all of the Langevin |
985 |
+ |
diffusion constants are within 6\% of the diffusion constant predicted |
986 |
+ |
from the fully solvated system. |
987 |
+ |
|
988 |
+ |
For rotational motion, Langevin integration (and the hydrodynamic tensor) |
989 |
+ |
yields rotational correlation times that are substantially shorter than those |
990 |
+ |
obtained from explicitly-solvated simulations. It is likely that this is due |
991 |
+ |
to the large size of the explicit solvent spheres, a feature that prevents |
992 |
+ |
the solvent from coming in contact with a substantial fraction of the surface |
993 |
+ |
area of the dumbbell. Therefore, the explicit solvent only provides drag |
994 |
+ |
over a substantially reduced surface area of this model, while the |
995 |
+ |
hydrodynamic theories utilize the entire surface area for estimating |
996 |
+ |
rotational diffusion. A sketch of the free volume available in the explicit |
997 |
+ |
solvent simulations is shown in figure \ref{fig:explanation}. |
998 |
+ |
|
999 |
+ |
|
1000 |
|
\begin{figure} |
1001 |
|
\centering |
1002 |
< |
\includegraphics[width=\linewidth]{vacf.eps} |
1003 |
< |
\caption[Plots of Velocity Auto-correlation Functions]{Velocity |
1004 |
< |
auto-correlation functions of NVE (explicit solvent) in blue and |
1005 |
< |
Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf} |
1002 |
> |
\includegraphics[width=6in]{explanation} |
1003 |
> |
\caption[Explanations of the differences between orientational |
1004 |
> |
correlation times for explicitly-solvated models and hydrodynamics |
1005 |
> |
predictions]{Explanations of the differences between orientational |
1006 |
> |
correlation times for explicitly-solvated models and hydrodynamic |
1007 |
> |
predictions. For the ellipsoids (upper figures), rotation of the |
1008 |
> |
principal axis can involve correlated collisions at both sides of the |
1009 |
> |
solute. In the rigid dumbbell model (lower figures), the large size |
1010 |
> |
of the explicit solvent spheres prevents them from coming in contact |
1011 |
> |
with a substantial fraction of the surface area of the dumbbell. |
1012 |
> |
Therefore, the explicit solvent only provides drag over a |
1013 |
> |
substantially reduced surface area of this model, where the |
1014 |
> |
hydrodynamic theories utilize the entire surface area for estimating |
1015 |
> |
rotational diffusion. |
1016 |
> |
} \label{fig:explanation} |
1017 |
|
\end{figure} |
1018 |
|
|
1019 |
+ |
|
1020 |
+ |
|
1021 |
+ |
\subsection{Composite banana-shaped molecules} |
1022 |
+ |
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids have |
1023 |
+ |
been used by Orlandi {\it et al.} to observe mesophases in |
1024 |
+ |
coarse-grained models for bent-core liquid crystalline |
1025 |
+ |
molecules.\cite{Orlandi:2006fk} We have used the same overlapping |
1026 |
+ |
ellipsoids as a way to test the behavior of our algorithm for a |
1027 |
+ |
structure of some interest to the materials science community, |
1028 |
+ |
although since we are interested in capturing only the hydrodynamic |
1029 |
+ |
behavior of this model, we have left out the dipolar interactions of |
1030 |
+ |
the original Orlandi model. |
1031 |
+ |
|
1032 |
+ |
A reference system composed of a single banana rigid body embedded in a |
1033 |
+ |
sea of 1929 solvent particles was created and run under standard |
1034 |
+ |
(microcanonical) molecular dynamics. The resulting viscosity of this |
1035 |
+ |
mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})). |
1036 |
+ |
To calculate the hydrodynamic properties of the banana rigid body model, |
1037 |
+ |
we created a rough shell (see Fig.~\ref{fig:roughShell}), in which |
1038 |
+ |
the banana is represented as a ``shell'' made of 3321 identical beads |
1039 |
+ |
(0.25 \AA\ in diameter) distributed on the surface. Applying the |
1040 |
+ |
procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
1041 |
+ |
identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 \AA), as |
1042 |
+ |
well as the resistance tensor, |
1043 |
+ |
\begin{equation*} |
1044 |
+ |
\Xi = |
1045 |
+ |
\left( {\begin{array}{*{20}c} |
1046 |
+ |
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
1047 |
+ |
0& 0.9270&-0.007063& 0.08585&0&0\\ |
1048 |
+ |
0&-0.007063&0.7494&0.2057&0&0\\ |
1049 |
+ |
0&0.0858&0.2057& 58.64& 0&0\\0.08585&0&0&0&48.30&3.219&\\0.2057&0&0&0&3.219&10.7373\\\end{array}} \right), |
1050 |
+ |
\end{equation*} |
1051 |
+ |
where the units for translational, translation-rotation coupling and |
1052 |
+ |
rotational tensors are (kcal fs / mol \AA$^2$), (kcal fs / mol \AA\ rad), |
1053 |
+ |
and (kcal fs / mol rad$^2$), respectively. |
1054 |
+ |
|
1055 |
+ |
The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor) |
1056 |
+ |
are essentially quantitative for translational diffusion of this model. |
1057 |
+ |
Orientational correlation times under the Langevin rigid-body integrator |
1058 |
+ |
are within 11\% of the values obtained from explicit solvent, but these |
1059 |
+ |
models also exhibit some solvent inaccessible surface area in the |
1060 |
+ |
explicitly-solvated case. |
1061 |
+ |
|
1062 |
+ |
\subsection{Composite sphero-ellipsoids} |
1063 |
+ |
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
1064 |
+ |
used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01} |
1065 |
+ |
|
1066 |
+ |
MORE DETAILS |
1067 |
+ |
|
1068 |
+ |
|
1069 |
+ |
\subsection{Summary} |
1070 |
+ |
According to our simulations, the langevin dynamics is a reliable |
1071 |
+ |
theory to apply to replace the explicit solvents, especially for the |
1072 |
+ |
translation properties. For large molecules, the rotation properties |
1073 |
+ |
are also mimiced reasonablly well. |
1074 |
+ |
|
1075 |
+ |
\begin{table*} |
1076 |
+ |
\begin{minipage}{\linewidth} |
1077 |
+ |
\begin{center} |
1078 |
+ |
\caption{Translational diffusion constants (D) for the model systems |
1079 |
+ |
calculated using microcanonical simulations (with explicit solvent), |
1080 |
+ |
theoretical predictions, and Langevin simulations (with implicit solvent). |
1081 |
+ |
Analytical solutions for the exactly-solved hydrodynamics models are |
1082 |
+ |
from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936} |
1083 |
+ |
(ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq} |
1084 |
+ |
(dumbbell). The other model systems have no known analytic solution. |
1085 |
+ |
All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= |
1086 |
+ |
$10^{-4}$ \AA$^2$ / fs). } |
1087 |
+ |
\begin{tabular}{lccccccc} |
1088 |
+ |
\hline |
1089 |
+ |
& \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\ |
1090 |
+ |
\cline{2-3} \cline{5-7} |
1091 |
+ |
model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\ |
1092 |
+ |
\hline |
1093 |
+ |
sphere & 0.261 & ? & & 2.59 & exact & 2.59 & 2.56 \\ |
1094 |
+ |
ellipsoid & 0.255 & 2.44 & & 2.34 & exact & 2.34 & 2.37 \\ |
1095 |
+ |
& 0.255 & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\ |
1096 |
+ |
dumbbell & 0.322 & ? & & 1.57 & bead model & 1.57 & 1.57 \\ |
1097 |
+ |
& 0.322 & ? & & 1.57 & rough shell & ? & ? \\ |
1098 |
+ |
banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\ |
1099 |
+ |
lipid & 0.349 & 0.96 & & & rough shell & 1.33 & 1.32 \\ |
1100 |
+ |
\end{tabular} |
1101 |
+ |
\label{tab:translation} |
1102 |
+ |
\end{center} |
1103 |
+ |
\end{minipage} |
1104 |
+ |
\end{table*} |
1105 |
+ |
|
1106 |
+ |
\begin{table*} |
1107 |
+ |
\begin{minipage}{\linewidth} |
1108 |
+ |
\begin{center} |
1109 |
+ |
\caption{Orientational relaxation times ($\tau$) for the model systems using |
1110 |
+ |
microcanonical simulation (with explicit solvent), theoretical |
1111 |
+ |
predictions, and Langevin simulations (with implicit solvent). All |
1112 |
+ |
relaxation times are for the rotational correlation function with |
1113 |
+ |
$\ell = 2$ and are reported in units of ps. The ellipsoidal model has |
1114 |
+ |
an exact solution for the orientational correlation time due to |
1115 |
+ |
Perrin, but the other model systems have no known analytic solution.} |
1116 |
+ |
\begin{tabular}{lccccccc} |
1117 |
+ |
\hline |
1118 |
+ |
& \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\ |
1119 |
+ |
\cline{2-3} \cline{5-7} |
1120 |
+ |
model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic & simulation \\ |
1121 |
+ |
\hline |
1122 |
+ |
sphere & 0.261 & & & 9.06 & exact & 9.06 & 9.11 \\ |
1123 |
+ |
ellipsoid & 0.255 & 46.7 & & 22.0 & exact & 22.0 & 22.2 \\ |
1124 |
+ |
& 0.255 & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\ |
1125 |
+ |
dumbbell & 0.322 & 14.0 & & & bead model & 52.3 & 52.8 \\ |
1126 |
+ |
& 0.322 & 14.0 & & & rough shell & ? & ? \\ |
1127 |
+ |
banana & 0.298 & 63.8 & & & rough shell & 70.9 & 70.9 \\ |
1128 |
+ |
lipid & 0.349 & 78.0 & & & rough shell & 76.9 & 77.9 \\ |
1129 |
+ |
\hline |
1130 |
+ |
\end{tabular} |
1131 |
+ |
\label{tab:rotation} |
1132 |
+ |
\end{center} |
1133 |
+ |
\end{minipage} |
1134 |
+ |
\end{table*} |
1135 |
+ |
|
1136 |
+ |
\section{Application: A rigid-body lipid bilayer} |
1137 |
+ |
|
1138 |
+ |
The Langevin dynamics integrator was applied to study the formation of |
1139 |
+ |
corrugated structures emerging from simulations of the coarse grained |
1140 |
+ |
lipid molecular models presented above. The initial configuration is |
1141 |
+ |
taken from our molecular dynamics studies on lipid bilayers with |
1142 |
+ |
lennard-Jones sphere solvents. The solvent molecules were excluded |
1143 |
+ |
from the system, and the experimental value for the viscosity of water |
1144 |
+ |
at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects |
1145 |
+ |
of the solvent. The absence of explicit solvent molecules and the |
1146 |
+ |
stability of the integrator allowed us to take timesteps of 50 fs. A |
1147 |
+ |
total simulation run time of 100 ns was sampled. |
1148 |
+ |
Fig. \ref{fig:bilayer} shows the configuration of the system after 100 |
1149 |
+ |
ns, and the ripple structure remains stable during the entire |
1150 |
+ |
trajectory. Compared with using explicit bead-model solvent |
1151 |
+ |
molecules, the efficiency of the simulation has increased by an order |
1152 |
+ |
of magnitude. |
1153 |
+ |
|
1154 |
|
\begin{figure} |
1155 |
|
\centering |
1156 |
< |
\includegraphics[width=\linewidth]{uacf.eps} |
1157 |
< |
\caption[Auto-correlation functions of the principle axis of the |
1158 |
< |
middle GB particle]{Auto-correlation functions of the principle axis |
1159 |
< |
of the middle GB particle of NVE (blue) and Langevin dynamics |
1160 |
< |
(red).} \label{langevin:uacf} |
1156 |
> |
\includegraphics[width=\linewidth]{bilayer} |
1157 |
> |
\caption[Snapshot of a bilayer of rigid-body models for lipids]{A |
1158 |
> |
snapshot of a bilayer composed of rigid-body models for lipid |
1159 |
> |
molecules evolving using the Langevin integrator described in this |
1160 |
> |
work.} \label{fig:bilayer} |
1161 |
|
\end{figure} |
1162 |
|
|
1163 |
|
\section{Conclusions} |
1164 |
|
|
1165 |
|
We have presented a new Langevin algorithm by incorporating the |
1166 |
|
hydrodynamics properties of arbitrary shaped molecules into an |
1167 |
< |
advanced symplectic integration scheme. The temperature control |
1168 |
< |
ability of this algorithm was demonstrated by a set of simulations |
1169 |
< |
with different viscosities. It was also shown to have significant |
1170 |
< |
advantage of producing rapid thermal equilibration over |
717 |
< |
Nos\'{e}-Hoover method. Further studies in systems involving banana |
718 |
< |
shaped molecules illustrated that the dynamic properties could be |
719 |
< |
preserved by using this new algorithm as an implicit solvent model. |
1167 |
> |
advanced symplectic integration scheme. Further studies in systems |
1168 |
> |
involving banana shaped molecules illustrated that the dynamic |
1169 |
> |
properties could be preserved by using this new algorithm as an |
1170 |
> |
implicit solvent model. |
1171 |
|
|
1172 |
|
|
1173 |
|
\section{Acknowledgments} |
1177 |
|
of Notre Dame. |
1178 |
|
\newpage |
1179 |
|
|
1180 |
< |
\bibliographystyle{jcp2} |
1180 |
> |
\bibliographystyle{jcp} |
1181 |
|
\bibliography{langevin} |
1182 |
|
|
1183 |
|
\end{document} |