--- trunk/tengDissertation/Langevin.tex 2006/06/11 02:23:00 2853 +++ trunk/tengDissertation/Langevin.tex 2006/06/23 21:33:52 2882 @@ -1,5 +1,5 @@ -\chapter{\label{chapt:methodology}LANGEVIN DYNAMICS for RIGID BODIES of ARBITRARY SHAPE} +\chapter{\label{chapt:langevin}LANGEVIN DYNAMICS FOR RIGID BODIES OF ARBITRARY SHAPE} \section{Introduction} @@ -98,16 +98,15 @@ sophisticated rigid body dynamics. estimation of friction tensor from hydrodynamics theory into the sophisticated rigid body dynamics. -\section{Method{\label{methodSec}}} +\section{Computational Methods{\label{methodSec}}} \subsection{\label{introSection:frictionTensor}Friction Tensor} -Theoretically, the friction kernel can be determined using velocity -autocorrelation function. However, this approach become impractical -when the system become more and more complicate. Instead, various -approaches based on hydrodynamics have been developed to calculate -the friction coefficients. The friction effect is isotropic in -Equation, $\zeta$ can be taken as a scalar. In general, friction -tensor $\Xi$ is a $6\times 6$ matrix given by +Theoretically, the friction kernel can be determined using the +velocity autocorrelation function. However, this approach become +impractical when the system become more and more complicate. +Instead, various approaches based on hydrodynamics have been +developed to calculate the friction coefficients. In general, +friction tensor $\Xi$ is a $6\times 6$ matrix given by \[ \Xi = \left( {\begin{array}{*{20}c} {\Xi _{}^{tt} } & {\Xi _{}^{rt} } \\ @@ -138,8 +137,9 @@ For a spherical particle, the translational and rotati \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}} -For a spherical particle, the translational and rotational friction -constant can be calculated from Stoke's law, +For a spherical particle with slip boundary conditions, the +translational and rotational friction constant can be calculated +from Stoke's law, \[ \Xi ^{tt} = \left( {\begin{array}{*{20}c} {6\pi \eta R} & 0 & 0 \\ @@ -161,9 +161,9 @@ exactly. In 1936, Perrin extended Stokes's law to gene Other non-spherical shape, such as cylinder and ellipsoid \textit{etc}, are widely used as reference for developing new hydrodynamics theory, because their properties can be calculated -exactly. In 1936, Perrin extended Stokes's law to general ellipsoid, -also called a triaxial ellipsoid, which is given in Cartesian -coordinates by\cite{Perrin1934, Perrin1936} +exactly. In 1936, Perrin extended Stokes's law to general +ellipsoids, also called a triaxial ellipsoid, which is given in +Cartesian coordinates by\cite{Perrin1934, Perrin1936} \[ \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2 }} = 1 @@ -198,20 +198,20 @@ and \end{array}. \] -\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}} +\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}} -Unlike spherical and other regular shaped molecules, there is not -analytical solution for friction tensor of any arbitrary shaped +Unlike spherical and other simply shaped molecules, there is no +analytical solution for the friction tensor for arbitrarily shaped rigid molecules. The ellipsoid of revolution model and general triaxial ellipsoid model have been used to approximate the hydrodynamic properties of rigid bodies. However, since the mapping -from all possible ellipsoidal space, $r$-space, to all possible -combination of rotational diffusion coefficients, $D$-space is not +from all possible ellipsoidal spaces, $r$-space, to all possible +combination of rotational diffusion coefficients, $D$-space, is not unique\cite{Wegener1979} as well as the intrinsic coupling between -translational and rotational motion of rigid body, general ellipsoid -is not always suitable for modeling arbitrarily shaped rigid -molecule. A number of studies have been devoted to determine the -friction tensor for irregularly shaped rigid bodies using more +translational and rotational motion of rigid body, general +ellipsoids are not always suitable for modeling arbitrarily shaped +rigid molecule. A number of studies have been devoted to determine +the friction tensor for irregularly shaped rigid bodies using more advanced method where the molecule of interest was modeled by combinations of spheres(beads)\cite{Carrasco1999} and the hydrodynamics properties of the molecule can be calculated using the @@ -273,8 +273,8 @@ where $\delta _{ij}$ is Kronecker delta function. Inve B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij} )T_{ij} \] -where $\delta _{ij}$ is Kronecker delta function. Inverting matrix -$B$, we obtain +where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$ +matrix, we obtain \[ C = B^{ - 1} = \left( {\begin{array}{*{20}c} @@ -306,20 +306,21 @@ resistance (reaction), at which the trace of rotationa The resistance tensor depends on the origin to which they refer. The proper location for applying friction force is the center of -resistance (reaction), at which the trace of rotational resistance -tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of -resistance is defined as an unique point of the rigid body at which -the translation-rotation coupling tensor are symmetric, +resistance (or center of reaction), at which the trace of rotational +resistance tensor, $ \Xi ^{rr}$ reaches a minimum value. +Mathematically, the center of resistance is defined as an unique +point of the rigid body at which the translation-rotation coupling +tensor are symmetric, \begin{equation} \Xi^{tr} = \left( {\Xi^{tr} } \right)^T \label{introEquation:definitionCR} \end{equation} -Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, +From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, we can easily find out that the translational resistance tensor is origin independent, while the rotational resistance tensor and translation-rotation coupling resistance tensor depend on the -origin. Given resistance tensor at an arbitrary origin $O$, and a -vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can +origin. Given the resistance tensor at an arbitrary origin $O$, and +a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can obtain the resistance tensor at $P$ by \begin{equation} \begin{array}{l} @@ -360,7 +361,7 @@ joining center of resistance $R$ and origin $O$. where $x_OR$, $y_OR$, $z_OR$ are the components of the vector joining center of resistance $R$ and origin $O$. -\subsection{Langevin dynamics for rigid particles of arbitrary shape} +\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} Consider a Langevin equation of motions in generalized coordinates \begin{equation} @@ -369,7 +370,7 @@ $V_i = V_i(v_i,\omega _i)$. The right side of Eq. \end{equation} where $M_i$ is a $6\times6$ generalized diagonal mass (include mass and moment of inertial) matrix and $V_i$ is a generalized velocity, -$V_i = V_i(v_i,\omega _i)$. The right side of Eq. +$V_i = V_i(v_i,\omega _i)$. The right side of Eq.~ (\ref{LDGeneralizedForm}) consists of three generalized forces in lab-fixed frame, systematic force $F_{s,i}$, dissipative force $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the @@ -403,7 +404,7 @@ with zero mean and variance \begin{equation} \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle = \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle = -2k_B T\Xi _R \delta (t - t'). +2k_B T\Xi _R \delta (t - t'). \label{randomForce} \end{equation} The equation of motion for $v_i$ can be written as @@ -515,6 +516,122 @@ be advanced to the same time value. + \frac{h}{2} {\bf \tau}^b(t + h) . \end{align*} -\section{Results and discussion} +\section{Results and Discussion} +The Langevin algorithm described in previous section has been +implemented in {\sc oopse}\cite{Meineke2005} and applied to the +studies of kinetics and thermodynamic properties in several systems. + +\subsection{Temperature Control} + +As shown in Eq.~\ref{randomForce}, random collisions associated with +the solvent's thermal motions is controlled by the external +temperature. The capability to maintain the temperature of the whole +system was usually used to measure the stability and efficiency of +the algorithm. In order to verify the stability of this new +algorithm, a series of simulations are performed on system +consisiting of 256 SSD water molecules with different viscosities. +Initial configuration for the simulations is taken from a 1ns NVT +simulation with a cubic box of 19.7166~\AA. All simulation are +carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns +with reference temperature at 300~K. Average temperature as a +function of $\eta$ is shown in Table \ref{langevin:viscosity} where +the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 - +1$ poise. The better temperature control at higher viscosity can be +explained by the finite size effect and relative slow relaxation +rate at lower viscosity regime. +\begin{table} +\caption{Average temperatures from Langevin dynamics simulations of +SSD water molecules with reference temperature at 300~K.} +\label{langevin:viscosity} +\begin{center} +\begin{tabular}{|l|l|l|} + \hline + $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\ + 1 & 300.47 & 10.99 \\ + 0.1 & 301.19 & 11.136 \\ + 0.01 & 303.04 & 11.796 \\ + \hline +\end{tabular} +\end{center} +\end{table} + +Another set of calculation were performed to study the efficiency of +temperature control using different temperature coupling schemes. +The starting configuration is cooled to 173~K and evolved using NVE, +NVT, and Langevin dynamic with time step of 2 fs. +Fig.~\ref{langevin:temperature} shows the heating curve obtained as +the systems reach equilibrium. The orange curve in +Fig.~\ref{langevin:temperature} represents the simulation using +Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps +which gives reasonable tight coupling, while the blue one from +Langevin dynamics with viscosity of 0.1 poise demonstrates a faster +scaling to the desire temperature. In extremely lower friction +regime (when $ \eta \approx 0$), Langevin dynamics becomes normal +NVE (see green curve in Fig.~\ref{langevin:temperature}) which loses +the temperature control ability. + +\begin{figure} +\centering +\includegraphics[width=\linewidth]{temperature.eps} +\caption[Plot of Temperature Fluctuation Versus Time]{Plot of +temperature fluctuation versus time.} \label{langevin:temperature} +\end{figure} + +\subsection{Langevin Dynamics of Banana Shaped Molecule} + +In order to verify that Langevin dynamics can mimic the dynamics of +the systems absent of explicit solvents, we carried out two sets of +simulations and compare their dynamic properties. + +\subsubsection{Simulations Contain One Banana Shaped Molecule} + +Fig.~\ref{langevin:oneBanana} shows a snapshot of the simulation +made of 256 pentane molecules and one banana shaped molecule at +273~K. It has an equivalent implicit solvent system containing only +one banana shaped molecule with viscosity of 0.289 center poise. To +calculate the hydrodynamic properties of the banana shaped molecule, +we create a rough shell model (see Fig.~\ref{langevin:roughShell}), +in which the banana shaped molecule is represented as a ``shell'' +made of 2266 small identical beads with size of 0.3 $\AA$ on the +surface. Applying the procedure described in +Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we +identified the center of resistance at $(0, 0.7482, -0.1988)$, as +well as the resistance tensor, +\[ +\left( {\begin{array}{*{20}c} +0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\ +3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\ +-6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\ +5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\ +0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\ +0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\ +\end{array}} \right). +\] + +\begin{figure} +\centering +\includegraphics[width=\linewidth]{roughShell.eps} +\caption[Rough shell model for banana shaped molecule]{Rough shell +model for banana shaped molecule.} \label{langevin:roughShell} +\end{figure} + +%\begin{figure} +%\centering +%\includegraphics[width=\linewidth]{oneBanana.eps} +%\caption[Snapshot from Simulation of One Banana Shaped Molecules and +%256 Pentane Molecules]{Snapshot from simulation of one Banana shaped +%molecules and 256 pentane molecules.} \label{langevin:oneBanana} +%\end{figure} + +\subsubsection{Simulations Contain Two Banana Shaped Molecules} + +\begin{figure} +\centering +\includegraphics[width=\linewidth]{twoBanana.eps} +\caption[Snapshot from Simulation of Two Banana Shaped Molecules and +256 Pentane Molecules]{Snapshot from simulation of two Banana shaped +molecules and 256 pentane molecules.} \label{langevin:twoBanana} +\end{figure} + \section{Conclusions}