| 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. |
| 173 |
– |
|
| 174 |
– |
\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 |
|
|
| 429 |
– |
\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} |
| 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 |
| 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 & & & \\ |
| 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 & & & \\ |
| 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 \\ |
| 644 |
|
\end{figure} |
| 645 |
|
|
| 646 |
|
\subsection{Simulation Methodology} |
| 634 |
– |
|
| 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 |
| 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. A switching function was applied to |
| 656 |
< |
all potentials to smoothly turn off the interactions between a range |
| 657 |
< |
of $22$ and $25$ \AA. The switching function was the standard (cubic) |
| 658 |
< |
function, |
| 655 |
> |
able to use a time step of 25 fs. |
| 656 |
> |
|
| 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 |
> |
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} |
| 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} \langle \left( |
| 721 |
< |
\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \rangle_{t_0}. |
| 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} \langle \left( |
| 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)-\langle P \rangle \right)dt' |
| 729 |
< |
\right)^2 \rangle_{t_0}. |
| 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} \langle P_{xz}(t_0) P_{xz}(t_0 |
| 735 |
< |
+ t) \rangle_{t_0} dt, |
| 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. |
| 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} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle, |
| 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 |
| 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) = \langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf |
| 823 |
< |
u}_{i}(0) \right) |
| 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 |
| 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 |
< |
qsolvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} |
| 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 |
| 865 |
|
$\Xi_{tt}$ assuming behavior intermediate between the two boundary |
| 866 |
|
conditions. |
| 867 |
|
|
| 868 |
< |
In these simulations, our spherical particles were structureless, so |
| 869 |
< |
there is no way to obtain rotational correlation times for this model |
| 870 |
< |
system. |
| 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 |
< |
\subsubsection{Ellipsoids} |
| 874 |
< |
For uniaxial ellipsoids ($a > b = c$) , Perrin's formulae for both |
| 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} |
| 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 modle ellipsoid. We |
| 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 |
| 914 |
|
exact treatment of the diffusion tensor as well as the rough-shell |
| 915 |
|
model for the ellipsoid. |
| 916 |
|
|
| 917 |
< |
The agreement with the translational diffusion constants from the |
| 918 |
< |
microcanonical simulations is quite good, although the rotational |
| 919 |
< |
correlation times are a factor of 2 shorter than the predictions of |
| 920 |
< |
the Perrin hydrodynamic model. |
| 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 |
< |
\subsubsection{Rigid dumbbells} |
| 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 |
| 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). |
| 961 |
> |
inversion of a 6 x 6 matrix). |
| 962 |
|
|
| 963 |
+ |
\begin{figure} |
| 964 |
+ |
\centering |
| 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 |
| 985 |
|
diffusion constants are within 6\% of the diffusion constant predicted |
| 986 |
|
from the fully solvated system. |
| 987 |
|
|
| 988 |
< |
For rotational motion, Langevin integration yields |
| 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 |
|
|
| 909 |
– |
\subsubsection{Ellipsoidal-composite banana-shaped molecules} |
| 999 |
|
|
| 1000 |
< |
Banana-shaped rigid bodies composed of composites of Gay-Berne |
| 1001 |
< |
ellipsoids have been used by Orlandi {\it et al.} to observe |
| 1002 |
< |
mesophases in coarse-grained models bent-core liquid crystalline |
| 1003 |
< |
molecules.\cite{Orlandi:2006fk} We have used the overlapping |
| 1000 |
> |
\begin{figure} |
| 1001 |
> |
\centering |
| 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 leave out the dipolar interactions of the |
| 1030 |
< |
original Orlandi model. |
| 1031 |
< |
|
| 1032 |
< |
\subsubsection{Composite sphero-ellipsoids} |
| 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{Temperature Control} |
| 929 |
< |
|
| 930 |
< |
As shown in Eq.~\ref{randomForce}, random collisions associated with |
| 931 |
< |
the solvent's thermal motions is controlled by the external |
| 932 |
< |
temperature. The capability to maintain the temperature of the whole |
| 933 |
< |
system was usually used to measure the stability and efficiency of |
| 934 |
< |
the algorithm. In order to verify the stability of this new |
| 935 |
< |
algorithm, a series of simulations are performed on system |
| 936 |
< |
consisiting of 256 SSD water molecules with different viscosities. |
| 937 |
< |
The initial configuration for the simulations is taken from a 1ns |
| 938 |
< |
NVT simulation with a cubic box of 19.7166~\AA. All simulation are |
| 939 |
< |
carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns |
| 940 |
< |
with reference temperature at 300~K. The average temperature as a |
| 941 |
< |
function of $\eta$ is shown in Table \ref{langevin:viscosity} where |
| 942 |
< |
the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 - |
| 943 |
< |
1$ poise. The better temperature control at higher viscosity can be |
| 944 |
< |
explained by the finite size effect and relative slow relaxation |
| 945 |
< |
rate at lower viscosity regime. |
| 946 |
< |
\begin{table} |
| 947 |
< |
\caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF |
| 948 |
< |
SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.} |
| 949 |
< |
\label{langevin:viscosity} |
| 950 |
< |
\begin{center} |
| 951 |
< |
\begin{tabular}{lll} |
| 952 |
< |
\hline |
| 953 |
< |
$\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\ |
| 954 |
< |
\hline |
| 955 |
< |
1 & 300.47 & 10.99 \\ |
| 956 |
< |
0.1 & 301.19 & 11.136 \\ |
| 957 |
< |
0.01 & 303.04 & 11.796 \\ |
| 958 |
< |
\hline |
| 959 |
< |
\end{tabular} |
| 960 |
< |
\end{center} |
| 961 |
< |
\end{table} |
| 962 |
< |
|
| 963 |
< |
Another set of calculations were performed to study the efficiency of |
| 964 |
< |
temperature control using different temperature coupling schemes. |
| 965 |
< |
The starting configuration is cooled to 173~K and evolved using NVE, |
| 966 |
< |
NVT, and Langevin dynamic with time step of 2 fs. |
| 967 |
< |
Fig.~\ref{langevin:temperature} shows the heating curve obtained as |
| 968 |
< |
the systems reach equilibrium. The orange curve in |
| 969 |
< |
Fig.~\ref{langevin:temperature} represents the simulation using |
| 970 |
< |
Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps |
| 971 |
< |
which gives reasonable tight coupling, while the blue one from |
| 972 |
< |
Langevin dynamics with viscosity of 0.1 poise demonstrates a faster |
| 973 |
< |
scaling to the desire temperature. When $ \eta = 0$, Langevin dynamics becomes normal |
| 974 |
< |
NVE (see orange curve in Fig.~\ref{langevin:temperature}) which |
| 975 |
< |
loses the temperature control ability. |
| 976 |
< |
|
| 977 |
< |
\begin{figure} |
| 978 |
< |
\centering |
| 979 |
< |
\includegraphics[width=\linewidth]{temperature} |
| 980 |
< |
\caption[Plot of Temperature Fluctuation Versus Time]{Plot of |
| 981 |
< |
temperature fluctuation versus time.} \label{langevin:temperature} |
| 982 |
< |
\end{figure} |
| 983 |
< |
|
| 984 |
< |
|
| 985 |
< |
The diffusion constants and rotation relaxation times for |
| 986 |
< |
different shaped molecules are shown in table \ref{tab:translation} |
| 987 |
< |
and \ref{tab:rotation} to compare to the results calculated from NVE |
| 988 |
< |
simulations. The theoretical values for sphere is calculated from the |
| 989 |
< |
Stokes-Einstein law, the theoretical values for ellipsoid is |
| 990 |
< |
calculated from Perrin's fomula, The exact method is |
| 991 |
< |
applied to the langevin dynamics simulations for sphere and ellipsoid, |
| 992 |
< |
the bead model is applied to the simulation for dumbbell molecule, and |
| 993 |
< |
the rough shell model is applied to ellipsoid, dumbbell, banana and |
| 994 |
< |
lipid molecules. The results from all the langevin dynamics |
| 995 |
< |
simulations, including exact, bead model and rough shell, match the |
| 996 |
< |
theoretical values perfectly for all different shaped molecules. This |
| 997 |
< |
indicates that our simulation package for langevin dynamics is working |
| 998 |
< |
well. The approxiate methods ( bead model and rough shell model) are |
| 999 |
< |
accurate enough for the current simulations. The goal of the langevin |
| 1000 |
< |
dynamics theory is to replace the explicit solvents by the friction |
| 1001 |
< |
forces. We compared the dynamic properties of different shaped |
| 1002 |
< |
molecules in langevin dynamics simulations with that in NVE |
| 1003 |
< |
simulations. The results are reasonable close. Overall, the |
| 1004 |
< |
translational diffusion constants calculated from langevin dynamics |
| 1005 |
< |
simulations are very close to the values from the NVE simulation. For |
| 1006 |
< |
sphere and lipid molecules, the diffusion constants are a little bit |
| 1007 |
< |
off from the NVE simulation results. One possible reason is that the |
| 1008 |
< |
calculation of the viscosity is very difficult to be accurate. Another |
| 1009 |
< |
possible reason is that although we save very frequently during the |
| 1010 |
< |
NVE simulations and run pretty long time simulations, there is only |
| 1011 |
< |
one solute molecule in the system which makes the calculation for the |
| 1012 |
< |
diffusion constant difficult. The sphere molecule behaves as a free |
| 1013 |
< |
rotor in the solvent, so there is no rotation relaxation time |
| 1014 |
< |
calculated from NVE simulations. The rotation relaxation time is not |
| 1015 |
< |
very close to the NVE simulations results. The banana and lipid |
| 1016 |
< |
molecules match the NVE simulations results pretty well. The mismatch |
| 1017 |
< |
between langevin dynamics and NVE simulation for ellipsoid is possibly |
| 1018 |
< |
caused by the slip boundary condition. For dumbbell, the mismatch is |
| 1019 |
< |
caused by the size of the solvent molecule is pretty large compared to |
| 1020 |
< |
dumbbell molecule in NVE simulations. |
| 1021 |
< |
|
| 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 |
| 1133 |
|
\end{minipage} |
| 1134 |
|
\end{table*} |
| 1135 |
|
|
| 1136 |
< |
Langevin dynamics simulations are applied to study the formation of |
| 1137 |
< |
the ripple phase of lipid membranes. The initial configuration is |
| 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 are excluded from |
| 1143 |
< |
the system, the experimental value of water viscosity is applied to |
| 1144 |
< |
mimic the heat bath. Fig. XXX is the snapshot of the stable |
| 1145 |
< |
configuration of the system, the ripple structure stayed stable after |
| 1146 |
< |
100 ns run. The efficiency of the simulation is increased by one order |
| 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 |
|
|
| 1098 |
– |
\subsection{Langevin Dynamics of Banana Shaped Molecules} |
| 1099 |
– |
|
| 1100 |
– |
In order to verify that Langevin dynamics can mimic the dynamics of |
| 1101 |
– |
the systems absent of explicit solvents, we carried out two sets of |
| 1102 |
– |
simulations and compare their dynamic properties. |
| 1103 |
– |
Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation |
| 1104 |
– |
made of 256 pentane molecules and two banana shaped molecules at |
| 1105 |
– |
273~K. It has an equivalent implicit solvent system containing only |
| 1106 |
– |
two banana shaped molecules with viscosity of 0.289 center poise. To |
| 1107 |
– |
calculate the hydrodynamic properties of the banana shaped molecule, |
| 1108 |
– |
we created a rough shell model (see Fig.~\ref{langevin:roughShell}), |
| 1109 |
– |
in which the banana shaped molecule is represented as a ``shell'' |
| 1110 |
– |
made of 2266 small identical beads with size of 0.3 \AA on the |
| 1111 |
– |
surface. Applying the procedure described in |
| 1112 |
– |
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
| 1113 |
– |
identified the center of resistance at (0 $\rm{\AA}$, 0.7482 $\rm{\AA}$, |
| 1114 |
– |
-0.1988 $\rm{\AA}$), as well as the resistance tensor, |
| 1115 |
– |
\[ |
| 1116 |
– |
\left( {\begin{array}{*{20}c} |
| 1117 |
– |
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
| 1118 |
– |
0& 0.9270&-0.007063& 0.08585&0&0\\ |
| 1119 |
– |
0&-0.007063&0.7494&0.2057&0&0\\ |
| 1120 |
– |
0&0.0858&0.2057& 58.64& 0&0\\ |
| 1121 |
– |
0.08585&0&0&0&48.30&3.219&\\ |
| 1122 |
– |
0.2057&0&0&0&3.219&10.7373\\ |
| 1123 |
– |
\end{array}} \right). |
| 1124 |
– |
\] |
| 1125 |
– |
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. |
| 1126 |
– |
Curves of the velocity auto-correlation functions in |
| 1127 |
– |
Fig.~\ref{langevin:vacf} were shown to match each other very well. |
| 1128 |
– |
However, because of the stochastic nature, simulation using Langevin |
| 1129 |
– |
dynamics was shown to decay slightly faster than MD. In order to |
| 1130 |
– |
study the rotational motion of the molecules, we also calculated the |
| 1131 |
– |
auto-correlation function of the principle axis of the second GB |
| 1132 |
– |
particle, $u$. The discrepancy shown in Fig.~\ref{langevin:uacf} was |
| 1133 |
– |
probably due to the reason that we used the experimental viscosity directly instead of calculating bulk viscosity from simulation. |
| 1134 |
– |
|
| 1154 |
|
\begin{figure} |
| 1155 |
|
\centering |
| 1156 |
< |
\includegraphics[width=\linewidth]{roughShell} |
| 1157 |
< |
\caption[Rough shell model for banana shaped molecule]{Rough shell |
| 1158 |
< |
model for banana shaped molecule.} \label{langevin:roughShell} |
| 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 |
|
|
| 1142 |
– |
\begin{figure} |
| 1143 |
– |
\centering |
| 1144 |
– |
\includegraphics[width=\linewidth]{twoBanana} |
| 1145 |
– |
\caption[Snapshot from Simulation of Two Banana Shaped Molecules and |
| 1146 |
– |
256 Pentane Molecules]{Snapshot from simulation of two Banana shaped |
| 1147 |
– |
molecules and 256 pentane molecules.} \label{langevin:twoBanana} |
| 1148 |
– |
\end{figure} |
| 1149 |
– |
|
| 1150 |
– |
\begin{figure} |
| 1151 |
– |
\centering |
| 1152 |
– |
\includegraphics[width=\linewidth]{vacf} |
| 1153 |
– |
\caption[Plots of Velocity Auto-correlation Functions]{Velocity |
| 1154 |
– |
auto-correlation functions of NVE (explicit solvent) in blue and |
| 1155 |
– |
Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf} |
| 1156 |
– |
\end{figure} |
| 1157 |
– |
|
| 1158 |
– |
\begin{figure} |
| 1159 |
– |
\centering |
| 1160 |
– |
\includegraphics[width=\linewidth]{uacf} |
| 1161 |
– |
\caption[Auto-correlation functions of the principle axis of the |
| 1162 |
– |
middle GB particle]{Auto-correlation functions of the principle axis |
| 1163 |
– |
of the middle GB particle of NVE (blue) and Langevin dynamics |
| 1164 |
– |
(red).} \label{langevin:uacf} |
| 1165 |
– |
\end{figure} |
| 1166 |
– |
|
| 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 |
| 1175 |
< |
Nos\'{e}-Hoover method. Further studies in systems involving banana |
| 1176 |
< |
shaped molecules illustrated that the dynamic properties could be |
| 1177 |
< |
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} |
