| 4 |
|
\usepackage{endfloat} |
| 5 |
|
\usepackage{amsmath,bm} |
| 6 |
|
\usepackage{amssymb} |
| 7 |
– |
\usepackage{epsf} |
| 7 |
|
\usepackage{times} |
| 8 |
|
\usepackage{mathptmx} |
| 9 |
|
\usepackage{setspace} |
| 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: |
| 27 |
|
gezelter@nd.edu} \\ |
| 631 |
|
|
| 632 |
|
\begin{figure} |
| 633 |
|
\centering |
| 634 |
< |
\includegraphics[width=\linewidth]{temperature.eps} |
| 634 |
> |
\includegraphics[width=\linewidth]{temperature.pdf} |
| 635 |
|
\caption[Plot of Temperature Fluctuation Versus Time]{Plot of |
| 636 |
|
temperature fluctuation versus time.} \label{langevin:temperature} |
| 637 |
|
\end{figure} |
| 638 |
+ |
|
| 639 |
+ |
\subsection{Langevin Dynamics simulation {\it vs} NVE simulations} |
| 640 |
|
|
| 641 |
+ |
To validate our langevin dynamics simulation. We performed several NVE |
| 642 |
+ |
simulations with explicit solvents for different shaped |
| 643 |
+ |
molecules. There are one solute molecule and 1929 solvent molecules in |
| 644 |
+ |
NVE simulation. The parameters are shown in table |
| 645 |
+ |
\ref{tab:parameters}. The force field between spheres is standard |
| 646 |
+ |
Lennard-Jones, and ellipsoids interact with other ellipsoids and |
| 647 |
+ |
spheres with generalized Gay-Berne potential. All simulations are |
| 648 |
+ |
carried out at 300 K and 1 Atm. The time step is 25 ns, and a |
| 649 |
+ |
switching function was applied to all potentials to smoothly turn off |
| 650 |
+ |
the interactions between a range of $22$ and $25$ \AA. The switching |
| 651 |
+ |
function was the standard (cubic) function, |
| 652 |
+ |
\begin{equation} |
| 653 |
+ |
s(r) = |
| 654 |
+ |
\begin{cases} |
| 655 |
+ |
1 & \text{if $r \le r_{\text{sw}}$},\\ |
| 656 |
+ |
\frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2} |
| 657 |
+ |
{(r_{\text{cut}} - r_{\text{sw}})^3} |
| 658 |
+ |
& \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\ |
| 659 |
+ |
0 & \text{if $r > r_{\text{cut}}$.} |
| 660 |
+ |
\end{cases} |
| 661 |
+ |
\label{eq:switchingFunc} |
| 662 |
+ |
\end{equation} |
| 663 |
+ |
We have computed translational diffusion constants for lipid molecules |
| 664 |
+ |
from the mean-square displacement, |
| 665 |
+ |
\begin{equation} |
| 666 |
+ |
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle, |
| 667 |
+ |
\end{equation} |
| 668 |
+ |
of the solute molecules. Translational diffusion constants for the |
| 669 |
+ |
different shaped molecules are shown in table |
| 670 |
+ |
\ref{tab:translation}. We have also computed orientational correlation |
| 671 |
+ |
times for different shaped molecules from fits of the second-order |
| 672 |
+ |
Legendre polynomial correlation function, |
| 673 |
+ |
\begin{equation} |
| 674 |
+ |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
| 675 |
+ |
\mu}_{i}(0) \right) |
| 676 |
+ |
\end{equation} |
| 677 |
+ |
the results are shown in table \ref{tab:rotation}. We used einstein |
| 678 |
+ |
format of the pressure correlation function, |
| 679 |
+ |
\begin{equation} |
| 680 |
+ |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
| 681 |
+ |
\mu}_{i}(0) \right) |
| 682 |
+ |
\end{equation} |
| 683 |
+ |
to estimate the viscosity of the systems from NVE simulations. The |
| 684 |
+ |
viscosity can also be calculated by Green-Kubo pressure correlaton |
| 685 |
+ |
function, |
| 686 |
+ |
\begin{equation} |
| 687 |
+ |
C_{\ell}(t) = \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf |
| 688 |
+ |
\mu}_{i}(0) \right) |
| 689 |
+ |
\end{equation} |
| 690 |
+ |
However, this method converges slowly, and the statistics are not good |
| 691 |
+ |
enough to give us a very accurate value. The langevin dynamics |
| 692 |
+ |
simulations for different shaped molecules are performed at the same |
| 693 |
+ |
conditions as the NVE simulations with viscosity estimated from NVE |
| 694 |
+ |
simulations. To get better statistics, 1024 non-interacting solute |
| 695 |
+ |
molecules are put into one simulation box for each langevin |
| 696 |
+ |
simulation, this is equal to 1024 simulations for single solute |
| 697 |
+ |
systems. The diffusion constants and rotation relaxation times for |
| 698 |
+ |
different shaped molecules are shown in table \ref{tab:translation} |
| 699 |
+ |
and \ref{tab:rotation} to compare to the results calculated from NVE |
| 700 |
+ |
simulations. The theoretical values for sphere is calculated from the |
| 701 |
+ |
Stokes-Einstein law, the theoretical values for ellipsoid is |
| 702 |
+ |
calculated from Perrin's fomula, the theoretical values for dumbbell |
| 703 |
+ |
is calculated from StinXX and Davis theory. The exact method is |
| 704 |
+ |
applied to the langevin dynamics simulations for sphere and ellipsoid, |
| 705 |
+ |
the bead model is applied to the simulation for dumbbell molecule, and |
| 706 |
+ |
the rough shell model is applied to ellipsoid, dumbbell, banana and |
| 707 |
+ |
lipid molecules. The results from all the langevin dynamics |
| 708 |
+ |
simulations, including exact, bead model and rough shell, match the |
| 709 |
+ |
theoretical values perfectly for all different shaped molecules. This |
| 710 |
+ |
indicates that our simulation package for langevin dynamics is working |
| 711 |
+ |
well. The approxiate methods ( bead model and rough shell model) are |
| 712 |
+ |
accurate enough for the current simulations. The goal of the langevin |
| 713 |
+ |
dynamics theory is to replace the explicit solvents by the friction |
| 714 |
+ |
forces. We compared the dynamic properties of different shaped |
| 715 |
+ |
molecules in langevin dynamics simulations with that in NVE |
| 716 |
+ |
simulations. The results are reasonable close. Overall, the |
| 717 |
+ |
translational diffusion constants calculated from langevin dynamics |
| 718 |
+ |
simulations are very close to the values from the NVE simulation. For |
| 719 |
+ |
sphere and lipid molecules, the diffusion constants are a little bit |
| 720 |
+ |
off from the NVE simulation results. One possible reason is that the |
| 721 |
+ |
calculation of the viscosity is very difficult to be accurate. Another |
| 722 |
+ |
possible reason is that although we save very frequently during the |
| 723 |
+ |
NVE simulations and run pretty long time simulations, there is only |
| 724 |
+ |
one solute molecule in the system which makes the calculation for the |
| 725 |
+ |
diffusion constant difficult. The sphere molecule behaves as a free |
| 726 |
+ |
rotor in the solvent, so there is no rotation relaxation time |
| 727 |
+ |
calculated from NVE simulations. The rotation relaxation time is not |
| 728 |
+ |
very close to the NVE simulations results. The banana and lipid |
| 729 |
+ |
molecules match the NVE simulations results pretty well. The mismatch |
| 730 |
+ |
between langevin dynamics and NVE simulation for ellipsoid is possibly |
| 731 |
+ |
caused by the slip boundary condition. For dumbbell, the mismatch is |
| 732 |
+ |
caused by the size of the solvent molecule is pretty large compared to |
| 733 |
+ |
dumbbell molecule in NVE simulations. |
| 734 |
+ |
|
| 735 |
+ |
According to our simulations, the langevin dynamics is a reliable |
| 736 |
+ |
theory to apply to replace the explicit solvents, especially for the |
| 737 |
+ |
translation properties. For large molecules, the rotation properties |
| 738 |
+ |
are also mimiced reasonablly well. |
| 739 |
+ |
|
| 740 |
+ |
\begin{table*} |
| 741 |
+ |
\begin{minipage}{\linewidth} |
| 742 |
+ |
\begin{center} |
| 743 |
+ |
\caption{} |
| 744 |
+ |
\begin{tabular}{llccccccc} |
| 745 |
+ |
\hline |
| 746 |
+ |
& & Sphere & Ellipsoid & Dumbbell(2 spheres) & Banana(3 ellpsoids) & |
| 747 |
+ |
Lipid(head) & lipid(tail) & Solvent \\ |
| 748 |
+ |
\hline |
| 749 |
+ |
$d$ (\AA) & & 6.5 & 4.6 & 6.5 & 4.2 & 6.5 & 4.6 & 4.7 \\ |
| 750 |
+ |
$l$ (\AA) & & $= d$ & 13.8 & $=d$ & 11.2 & $=d$ & 13.8 & 4.7 \\ |
| 751 |
+ |
$\epsilon^s$ (kcal/mol) & & 0.8 & 0.8 & 0.8 & 0.8 & 0.185 & 0.8 & 0.8 \\ |
| 752 |
+ |
$\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 & 1 & 0.2 & 1 & 0.2 & 1 \\ |
| 753 |
+ |
$m$ (amu) & & 190 & 200 & 190 & 240 & 196 & 760 & 72.06 \\ |
| 754 |
+ |
%$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\ |
| 755 |
+ |
%\multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\ |
| 756 |
+ |
%\multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\ |
| 757 |
+ |
%\multicolumn{2}c{$I_{zz}$} & 0 & 9000 & N/A \\ |
| 758 |
+ |
%$\mu$ (Debye) & & varied & 0 & 0 \\ |
| 759 |
+ |
\end{tabular} |
| 760 |
+ |
\label{tab:parameters} |
| 761 |
+ |
\end{center} |
| 762 |
+ |
\end{minipage} |
| 763 |
+ |
\end{table*} |
| 764 |
+ |
|
| 765 |
+ |
\begin{table*} |
| 766 |
+ |
\begin{minipage}{\linewidth} |
| 767 |
+ |
\begin{center} |
| 768 |
+ |
\caption{} |
| 769 |
+ |
\begin{tabular}{lccccc} |
| 770 |
+ |
\hline |
| 771 |
+ |
& & & & &Translation \\ |
| 772 |
+ |
\hline |
| 773 |
+ |
& NVE & & Theoretical & Langevin & \\ |
| 774 |
+ |
\hline |
| 775 |
+ |
& $\eta$ & D & D & method & D \\ |
| 776 |
+ |
\hline |
| 777 |
+ |
sphere & 3.480159e-03 & 1.643135e-04 & 1.942779e-04 & exact & 1.982283e-04 \\ |
| 778 |
+ |
ellipsoid & 2.551262e-03 & 2.437492e-04 & 2.335756e-04 & exact & 2.374905e-04 \\ |
| 779 |
+ |
& 2.551262e-03 & 2.437492e-04 & 2.335756e-04 & rough shell & 2.284088e-04 \\ |
| 780 |
+ |
dumbell & 2.41276e-03 & 2.129432e-04 & 2.090239e-04 & bead model & 2.148098e-04 \\ |
| 781 |
+ |
& 2.41276e-03 & 2.129432e-04 & 2.090239e-04 & rough shell & 2.013219e-04 \\ |
| 782 |
+ |
banana & 2.9846e-03 & 1.527819e-04 & & rough shell & 1.54807e-04 \\ |
| 783 |
+ |
lipid & 3.488661e-03 & 0.9562979e-04 & & rough shell & 1.320987e-04 \\ |
| 784 |
+ |
\end{tabular} |
| 785 |
+ |
\label{tab:translation} |
| 786 |
+ |
\end{center} |
| 787 |
+ |
\end{minipage} |
| 788 |
+ |
\end{table*} |
| 789 |
+ |
|
| 790 |
+ |
\begin{table*} |
| 791 |
+ |
\begin{minipage}{\linewidth} |
| 792 |
+ |
\begin{center} |
| 793 |
+ |
\caption{} |
| 794 |
+ |
\begin{tabular}{lccccc} |
| 795 |
+ |
\hline |
| 796 |
+ |
& & & & &Rotation \\ |
| 797 |
+ |
\hline |
| 798 |
+ |
& NVE & & Theoretical & Langevin & \\ |
| 799 |
+ |
\hline |
| 800 |
+ |
& $\eta$ & $\tau_0$ & $\tau_0$ & method & $\tau_0$ \\ |
| 801 |
+ |
\hline |
| 802 |
+ |
sphere & 3.480159e-03 & & 1.208178e+04 & exact & 1.20628e+04 \\ |
| 803 |
+ |
ellipsoid & 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & exact & 2.21507e+04 \\ |
| 804 |
+ |
& 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & rough shell & 2.21714e+04 \\ |
| 805 |
+ |
dumbell & 2.41276e-03 & 1.42974e+04 & & bead model & 7.12435e+04 \\ |
| 806 |
+ |
& 2.41276e-03 & 1.42974e+04 & & rough shell & 7.04765e+04 \\ |
| 807 |
+ |
banana & 2.9846e-03 & 6.38323e+04 & & rough shell & 7.0945e+04 \\ |
| 808 |
+ |
lipid & 3.488661e-03 & 7.79595e+04 & & rough shell & 7.78886e+04 \\ |
| 809 |
+ |
\end{tabular} |
| 810 |
+ |
\label{tab:rotation} |
| 811 |
+ |
\end{center} |
| 812 |
+ |
\end{minipage} |
| 813 |
+ |
\end{table*} |
| 814 |
+ |
|
| 815 |
+ |
Langevin dynamics simulations are applied to study the formation of |
| 816 |
+ |
the ripple phase of lipid membranes. The initial configuration is |
| 817 |
+ |
taken from our molecular dynamics studies on lipid bilayers with |
| 818 |
+ |
lennard-Jones sphere solvents. The solvent molecules are excluded from |
| 819 |
+ |
the system, the experimental value of water viscosity is applied to |
| 820 |
+ |
mimic the heat bath. Fig. XXX is the snapshot of the stable |
| 821 |
+ |
configuration of the system, the ripple structure stayed stable after |
| 822 |
+ |
100 ns run. The efficiency of the simulation is increased by one order |
| 823 |
+ |
of magnitude. |
| 824 |
+ |
|
| 825 |
|
\subsection{Langevin Dynamics of Banana Shaped Molecules} |
| 826 |
|
|
| 827 |
|
In order to verify that Langevin dynamics can mimic the dynamics of |
| 861 |
|
|
| 862 |
|
\begin{figure} |
| 863 |
|
\centering |
| 864 |
< |
\includegraphics[width=\linewidth]{roughShell.eps} |
| 864 |
> |
\includegraphics[width=\linewidth]{roughShell.pdf} |
| 865 |
|
\caption[Rough shell model for banana shaped molecule]{Rough shell |
| 866 |
|
model for banana shaped molecule.} \label{langevin:roughShell} |
| 867 |
|
\end{figure} |
| 868 |
|
|
| 869 |
|
\begin{figure} |
| 870 |
|
\centering |
| 871 |
< |
\includegraphics[width=\linewidth]{twoBanana.eps} |
| 871 |
> |
\includegraphics[width=\linewidth]{twoBanana.pdf} |
| 872 |
|
\caption[Snapshot from Simulation of Two Banana Shaped Molecules and |
| 873 |
|
256 Pentane Molecules]{Snapshot from simulation of two Banana shaped |
| 874 |
|
molecules and 256 pentane molecules.} \label{langevin:twoBanana} |
| 876 |
|
|
| 877 |
|
\begin{figure} |
| 878 |
|
\centering |
| 879 |
< |
\includegraphics[width=\linewidth]{vacf.eps} |
| 879 |
> |
\includegraphics[width=\linewidth]{vacf.pdf} |
| 880 |
|
\caption[Plots of Velocity Auto-correlation Functions]{Velocity |
| 881 |
|
auto-correlation functions of NVE (explicit solvent) in blue and |
| 882 |
|
Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf} |
| 884 |
|
|
| 885 |
|
\begin{figure} |
| 886 |
|
\centering |
| 887 |
< |
\includegraphics[width=\linewidth]{uacf.eps} |
| 887 |
> |
\includegraphics[width=\linewidth]{uacf.pdf} |
| 888 |
|
\caption[Auto-correlation functions of the principle axis of the |
| 889 |
|
middle GB particle]{Auto-correlation functions of the principle axis |
| 890 |
|
of the middle GB particle of NVE (blue) and Langevin dynamics |
