| 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 |
| 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 |
| 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 |
< |
\begin{figure} |
| 870 |
< |
\centering |
| 871 |
< |
\includegraphics[width=6in]{explanation} |
| 872 |
< |
\caption[Explanations of the differences between orientational correlation times for explicitly-solvated models and hydrodynamics predictions]{Explanations of the differences between orientational correlation times for explicitly-solvated models and hydrodynamic predictions. For the ellipsoids (upper figures), rotation of the principal axis can involve correlated collisions at both sides of the solute. In the rigid dumbbell model (lower figures), the large size of the explicit solvent spheres prevents them from coming in contact with a substantial fraction of the surface area of the dumbbell. Therefore, the explicit solvent only provides drag over a substantially reduced surface area of this model, where the hydrodynamic theories utilize the entire surface area for estimating rotational diffusion. |
| 873 |
< |
} \label{fig:explanation} |
| 874 |
< |
\end{figure} |
| 875 |
< |
|
| 876 |
< |
\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 |
| 962 |
|
|
| 963 |
|
\begin{figure} |
| 964 |
|
\centering |
| 965 |
< |
\includegraphics[width=3in]{RoughShell} |
| 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 |
| 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 |
< |
\subsubsection{Ellipsoidal-composite banana-shaped molecules} |
| 1000 |
< |
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids |
| 1001 |
< |
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 same overlapping |
| 999 |
> |
|
| 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 have left out the dipolar interactions of the |
| 1030 |
< |
original Orlandi model. |
| 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{langevin:roughShell}), in which |
| 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 |
| 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 at (0 $\rm{\AA}$, 0.81 $\rm{\AA}$, 0 $\rm{\AA}$), as well as the resistance tensor, |
| 1042 |
< |
\[ |
| 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 |
< |
\] |
| 1051 |
< |
where the units for translational, translation-rotation coupling and rotational |
| 1052 |
< |
tensors are $\frac{kcal \cdot fs}{mol \cdot \rm{\AA}^2}$, $\frac{kcal \cdot fs}{ |
| 1053 |
< |
mol \cdot \rm{\AA} \cdot rad}$ and $\frac{kcal \cdot fs}{mol \cdot rad^2}$ respe |
| 977 |
< |
ctively. |
| 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. |
| 1059 |
|
models also exhibit some solvent inaccessible surface area in the |
| 1060 |
|
explicitly-solvated case. |
| 1061 |
|
|
| 1062 |
< |
\subsubsection{Composite sphero-ellipsoids} |
| 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 |
< |
The diffusion constants and rotation relaxation times for |
| 991 |
< |
different shaped molecules are shown in table \ref{tab:translation} |
| 992 |
< |
and \ref{tab:rotation} to compare to the results calculated from NVE |
| 993 |
< |
simulations. The theoretical values for sphere is calculated from the |
| 994 |
< |
Stokes-Einstein law, the theoretical values for ellipsoid is |
| 995 |
< |
calculated from Perrin's fomula, The exact method is |
| 996 |
< |
applied to the langevin dynamics simulations for sphere and ellipsoid, |
| 997 |
< |
the bead model is applied to the simulation for dumbbell molecule, and |
| 998 |
< |
the rough shell model is applied to ellipsoid, dumbbell, banana and |
| 999 |
< |
lipid molecules. The results from all the langevin dynamics |
| 1000 |
< |
simulations, including exact, bead model and rough shell, match the |
| 1001 |
< |
theoretical values perfectly for all different shaped molecules. This |
| 1002 |
< |
indicates that our simulation package for langevin dynamics is working |
| 1003 |
< |
well. The approxiate methods ( bead model and rough shell model) are |
| 1004 |
< |
accurate enough for the current simulations. The goal of the langevin |
| 1005 |
< |
dynamics theory is to replace the explicit solvents by the friction |
| 1006 |
< |
forces. We compared the dynamic properties of different shaped |
| 1007 |
< |
molecules in langevin dynamics simulations with that in NVE |
| 1008 |
< |
simulations. The results are reasonable close. Overall, the |
| 1009 |
< |
translational diffusion constants calculated from langevin dynamics |
| 1010 |
< |
simulations are very close to the values from the NVE simulation. For |
| 1011 |
< |
sphere and lipid molecules, the diffusion constants are a little bit |
| 1012 |
< |
off from the NVE simulation results. One possible reason is that the |
| 1013 |
< |
calculation of the viscosity is very difficult to be accurate. Another |
| 1014 |
< |
possible reason is that although we save very frequently during the |
| 1015 |
< |
NVE simulations and run pretty long time simulations, there is only |
| 1016 |
< |
one solute molecule in the system which makes the calculation for the |
| 1017 |
< |
diffusion constant difficult. The sphere molecule behaves as a free |
| 1018 |
< |
rotor in the solvent, so there is no rotation relaxation time |
| 1019 |
< |
calculated from NVE simulations. The rotation relaxation time is not |
| 1020 |
< |
very close to the NVE simulations results. The banana and lipid |
| 1021 |
< |
molecules match the NVE simulations results pretty well. The mismatch |
| 1022 |
< |
between langevin dynamics and NVE simulation for ellipsoid is possibly |
| 1023 |
< |
caused by the slip boundary condition. For dumbbell, the mismatch is |
| 1024 |
< |
caused by the size of the solvent molecule is pretty large compared to |
| 1025 |
< |
dumbbell molecule in NVE simulations. |
| 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 |
| 1090 |
|
\cline{2-3} \cline{5-7} |
| 1091 |
|
model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\ |
| 1092 |
|
\hline |
| 1093 |
< |
sphere & 0.348 & 1.64 & & 1.94 & exact & 1.94 & 1.98 \\ |
| 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.241 & 2.13 & & 2.09 & bead model & 2.10 & 2.15 \\ |
| 1097 |
< |
& 0.241 & 2.13 & & 2.09 & rough shell & 2.03 & 2.01 \\ |
| 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} |
| 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.241 & 14.3 & & & bead model & 39.2 & 71.2 \\ |
| 1126 |
< |
& 0.241 & 14.3 & & & rough shell & 32.6 & 70.5 \\ |
| 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 |
| 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 |
|
|
| 1154 |
+ |
\begin{figure} |
| 1155 |
+ |
\centering |
| 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 |
