| 608 |
|
& & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\ |
| 609 |
|
& & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ (kcal/mol) & $\epsilon^r$ & |
| 610 |
|
$m$ (amu) & $I_{xx}$ & $I_{yy}$ & $I_{zz}$ \\ \hline |
| 611 |
< |
Sphere & & 6.5 & $= d$ & 0.8 & 1 & 190 & & & \\ |
| 611 |
> |
Sphere & & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
| 612 |
|
Ellipsoid & & 4.6 & 13.8 & 0.8 & 0.2 & 200 & 2105 & 2105 & 421 \\ |
| 613 |
< |
Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1 & 190 & & & \\ |
| 613 |
> |
Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75 & 802.75 \\ |
| 614 |
|
Banana &(3 identical ellipsoids)& 4.2 & 11.2 & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\ |
| 615 |
|
Lipid: & Spherical Head & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\ |
| 616 |
|
& Ellipsoidal Tail & 4.6 & 13.8 & 0.8 & 0.2 & 760 & 45000 & 45000 & 9000 \\ |
| 854 |
|
exact treatment of the diffusion tensor as well as the rough-shell |
| 855 |
|
model for the ellipsoid. |
| 856 |
|
|
| 857 |
< |
The agreement with the translational diffusion constants from the |
| 858 |
< |
microcanonical simulations is quite good, although the rotational |
| 859 |
< |
correlation times are a factor of 2 shorter than the predictions of |
| 860 |
< |
the Perrin hydrodynamic model. |
| 857 |
> |
The translational diffusion constants from the microcanonical simulations |
| 858 |
> |
agree well with the predictions of the Perrin model, although the rotational |
| 859 |
> |
correlation times are a factor of 2 shorter than expected from hydrodynamic |
| 860 |
> |
theory. One explanation for the slower rotation |
| 861 |
> |
of explicitly-solvated ellipsoids is the possibility that solute-solvent |
| 862 |
> |
collisions happen at both ends of the solute whenever the principal |
| 863 |
> |
axis of the ellipsoid is turning. In the upper portion of figure |
| 864 |
> |
\ref{fig:explanation} we sketch a physical picture of this explanation. |
| 865 |
> |
Since our Langevin integrator is providing nearly quantitative agreement with |
| 866 |
> |
the Perrin model, it also predicts orientational diffusion for ellipsoids that |
| 867 |
> |
exceed explicitly solvated correlation times by a factor of two. |
| 868 |
|
|
| 869 |
+ |
\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} |
| 877 |
|
Perhaps the only {\it composite} rigid body for which analytic |
| 878 |
|
expressions for the hydrodynamic tensor are available is the |
| 905 |
|
hydrodynamic tensor for a rough shell model can be quite expensive (in |
| 906 |
|
this case it requires inverting a 10104 x 10104 matrix), while the |
| 907 |
|
bead model is typically easy to compute (in this case requiring |
| 908 |
< |
inversion of a 6 x 6 matrix). |
| 908 |
> |
inversion of a 6 x 6 matrix). |
| 909 |
> |
|
| 910 |
> |
\begin{figure} |
| 911 |
> |
\centering |
| 912 |
> |
\includegraphics[width=3in]{RoughShell} |
| 913 |
> |
\caption[Model rigid bodies and their rough shell approximations]{The |
| 914 |
> |
model rigid bodies (left column) used to test this algorithm and their |
| 915 |
> |
rough-shell approximations (right-column) that were used to compute |
| 916 |
> |
the hydrodynamic tensors. The top two models (ellipsoid and dumbbell) |
| 917 |
> |
have analytic solutions and were used to test the rough shell |
| 918 |
> |
approximation. The lower two models (banana and lipid) were compared |
| 919 |
> |
with explicitly-solvated molecular dynamics simulations. } |
| 920 |
> |
\label{fig:roughShell} |
| 921 |
> |
\end{figure} |
| 922 |
|
|
| 923 |
+ |
|
| 924 |
|
Once the hydrodynamic tensor has been computed, there is no additional |
| 925 |
|
penalty for carrying out a Langevin simulation with either of the two |
| 926 |
|
different hydrodynamics models. Our naive expectation is that since |
| 932 |
|
diffusion constants are within 6\% of the diffusion constant predicted |
| 933 |
|
from the fully solvated system. |
| 934 |
|
|
| 935 |
< |
For rotational motion, Langevin integration yields |
| 935 |
> |
For rotational motion, Langevin integration (and the hydrodynamic tensor) |
| 936 |
> |
yields rotational correlation times that are substantially shorter than those |
| 937 |
> |
obtained from explicitly-solvated simulations. It is likely that this is due |
| 938 |
> |
to the large size of the explicit solvent spheres, a feature that prevents |
| 939 |
> |
the solvent from coming in contact with a substantial fraction of the surface |
| 940 |
> |
area of the dumbbell. Therefore, the explicit solvent only provides drag |
| 941 |
> |
over a substantially reduced surface area of this model, while the |
| 942 |
> |
hydrodynamic theories utilize the entire surface area for estimating |
| 943 |
> |
rotational diffusion. A sketch of the free volume available in the explicit |
| 944 |
> |
solvent simulations is shown in figure \ref{fig:explanation}. |
| 945 |
|
|
| 946 |
|
\subsubsection{Ellipsoidal-composite banana-shaped molecules} |
| 947 |
< |
|
| 948 |
< |
Banana-shaped rigid bodies composed of composites of Gay-Berne |
| 912 |
< |
ellipsoids have been used by Orlandi {\it et al.} to observe |
| 947 |
> |
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids |
| 948 |
> |
have been used by Orlandi {\it et al.} to observe |
| 949 |
|
mesophases in coarse-grained models bent-core liquid crystalline |
| 950 |
< |
molecules.\cite{Orlandi:2006fk} We have used the overlapping |
| 950 |
> |
molecules.\cite{Orlandi:2006fk} We have used the same overlapping |
| 951 |
|
ellipsoids as a way to test the behavior of our algorithm for a |
| 952 |
|
structure of some interest to the materials science community, |
| 953 |
|
although since we are interested in capturing only the hydrodynamic |
| 954 |
< |
behavior of this model, we leave out the dipolar interactions of the |
| 954 |
> |
behavior of this model, we have left out the dipolar interactions of the |
| 955 |
|
original Orlandi model. |
| 956 |
+ |
|
| 957 |
+ |
A reference system composed of a single banana rigid body embedded in a |
| 958 |
+ |
sea of 1929 solvent particles was created and run under standard |
| 959 |
+ |
(microcanonical) molecular dynamics. The resulting viscosity of this |
| 960 |
+ |
mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})). |
| 961 |
+ |
To calculate the hydrodynamic properties of the banana rigid body model, |
| 962 |
+ |
we created a rough shell (see Fig.~\ref{langevin:roughShell}), in which |
| 963 |
+ |
the banana is represented as a ``shell'' made of 3321 identical beads |
| 964 |
+ |
(0.25 \AA in diameter) distributed on the surface. Applying the |
| 965 |
+ |
procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we |
| 966 |
+ |
identified the center of resistance at (0 $\rm{\AA}$, 0.81 $\rm{\AA}$, 0 $\rm{\AA}$), as well as the resistance tensor, |
| 967 |
+ |
\[ |
| 968 |
+ |
\left( {\begin{array}{*{20}c} |
| 969 |
+ |
0.9261 & 0 & 0&0&0.08585&0.2057\\ |
| 970 |
+ |
0& 0.9270&-0.007063& 0.08585&0&0\\ |
| 971 |
+ |
0&-0.007063&0.7494&0.2057&0&0\\ |
| 972 |
+ |
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). |
| 973 |
+ |
\] |
| 974 |
+ |
where the units for translational, translation-rotation coupling and rotational |
| 975 |
+ |
tensors are $\frac{kcal \cdot fs}{mol \cdot \rm{\AA}^2}$, $\frac{kcal \cdot fs}{ |
| 976 |
+ |
mol \cdot \rm{\AA} \cdot rad}$ and $\frac{kcal \cdot fs}{mol \cdot rad^2}$ respe |
| 977 |
+ |
ctively. |
| 978 |
|
|
| 979 |
< |
\subsubsection{Composite sphero-ellipsoids} |
| 979 |
> |
The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor) |
| 980 |
> |
are essentially quantitative for translational diffusion of this model. |
| 981 |
> |
Orientational correlation times under the Langevin rigid-body integrator |
| 982 |
> |
are within 11\% of the values obtained from explicit solvent, but these |
| 983 |
> |
models also exhibit some solvent inaccessible surface area in the |
| 984 |
> |
explicitly-solvated case. |
| 985 |
|
|
| 986 |
+ |
\subsubsection{Composite sphero-ellipsoids} |
| 987 |
|
Spherical heads perched on the ends of Gay-Berne ellipsoids have been |
| 988 |
|
used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01} |
| 925 |
– |
|
| 926 |
– |
|
| 927 |
– |
|
| 928 |
– |
\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} |
| 989 |
|
|
| 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 |
– |
|
| 990 |
|
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 |
| 1047 |
|
\cline{2-3} \cline{5-7} |
| 1048 |
|
model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\ |
| 1049 |
|
\hline |
| 1050 |
< |
sphere & 0.261 & ? & & 2.59 & exact & 2.59 & 2.56 \\ |
| 1050 |
> |
sphere & 0.348 & 1.64 & & 1.94 & exact & 1.94 & 1.98 \\ |
| 1051 |
|
ellipsoid & 0.255 & 2.44 & & 2.34 & exact & 2.34 & 2.37 \\ |
| 1052 |
|
& 0.255 & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\ |
| 1053 |
< |
dumbbell & 0.322 & ? & & 1.57 & bead model & 1.57 & 1.57 \\ |
| 1054 |
< |
& 0.322 & ? & & 1.57 & rough shell & ? & ? \\ |
| 1053 |
> |
dumbbell & 0.241 & 2.13 & & 2.09 & bead model & 2.10 & 2.15 \\ |
| 1054 |
> |
& 0.241 & 2.13 & & 2.09 & rough shell & 2.03 & 2.01 \\ |
| 1055 |
|
banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\ |
| 1056 |
|
lipid & 0.349 & 0.96 & & & rough shell & 1.33 & 1.32 \\ |
| 1057 |
|
\end{tabular} |
| 1076 |
|
\cline{2-3} \cline{5-7} |
| 1077 |
|
model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic & simulation \\ |
| 1078 |
|
\hline |
| 1074 |
– |
sphere & 0.261 & & & 9.06 & exact & 9.06 & 9.11 \\ |
| 1079 |
|
ellipsoid & 0.255 & 46.7 & & 22.0 & exact & 22.0 & 22.2 \\ |
| 1080 |
|
& 0.255 & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\ |
| 1081 |
< |
dumbbell & 0.322 & 14.0 & & & bead model & 52.3 & 52.8 \\ |
| 1082 |
< |
& 0.322 & 14.0 & & & rough shell & ? & ? \\ |
| 1081 |
> |
dumbbell & 0.241 & 14.3 & & & bead model & 39.2 & 71.2 \\ |
| 1082 |
> |
& 0.241 & 14.3 & & & rough shell & 32.6 & 70.5 \\ |
| 1083 |
|
banana & 0.298 & 63.8 & & & rough shell & 70.9 & 70.9 \\ |
| 1084 |
|
lipid & 0.349 & 78.0 & & & rough shell & 76.9 & 77.9 \\ |
| 1085 |
|
\hline |
| 1099 |
|
100 ns run. The efficiency of the simulation is increased by one order |
| 1100 |
|
of magnitude. |
| 1101 |
|
|
| 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 |
– |
|
| 1135 |
– |
\begin{figure} |
| 1136 |
– |
\centering |
| 1137 |
– |
\includegraphics[width=\linewidth]{roughShell} |
| 1138 |
– |
\caption[Rough shell model for banana shaped molecule]{Rough shell |
| 1139 |
– |
model for banana shaped molecule.} \label{langevin:roughShell} |
| 1140 |
– |
\end{figure} |
| 1141 |
– |
|
| 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 |
– |
|
| 1102 |
|
\section{Conclusions} |
| 1103 |
|
|
| 1104 |
|
We have presented a new Langevin algorithm by incorporating the |
| 1105 |
|
hydrodynamics properties of arbitrary shaped molecules into an |
| 1106 |
< |
advanced symplectic integration scheme. The temperature control |
| 1107 |
< |
ability of this algorithm was demonstrated by a set of simulations |
| 1108 |
< |
with different viscosities. It was also shown to have significant |
| 1109 |
< |
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. |
| 1106 |
> |
advanced symplectic integration scheme. Further studies in systems |
| 1107 |
> |
involving banana shaped molecules illustrated that the dynamic |
| 1108 |
> |
properties could be preserved by using this new algorithm as an |
| 1109 |
> |
implicit solvent model. |
| 1110 |
|
|
| 1111 |
|
|
| 1112 |
|
\section{Acknowledgments} |
