[ViewVC] Diff of: group/trunk/tengDissertation/Introduction.tex

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2718 by tim, Tue Apr 18 04:11:56 2006 UTC vs.
Revision 2725 by tim, Fri Apr 21 05:45:14 2006 UTC

+q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
+q(\Delta t)} \right]. %
-<
+ \label{introEquation:positionVerlet1}
->
+ \label{introEquation:positionVerlet2}
+\end{align}
+\subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
+\section{\label{introSection:molecularDynamics}Molecular Dynamics}
-<
+As a special discipline of molecular modeling, Molecular dynamics
-<
+has proven to be a powerful tool for studying the functions of
-<
+biological systems, providing structural, thermodynamic and
-<
+dynamical information.
-<
-<
+\subsection{\label{introSec:mdInit}Initialization}
->
+As one of the principal tools of molecular modeling, Molecular
->
+dynamics has proven to be a powerful tool for studying the functions
->
+of biological systems, providing structural, thermodynamic and
->
+dynamical information. The basic idea of molecular dynamics is that
->
+macroscopic properties are related to microscopic behavior and
->
+microscopic behavior can be calculated from the trajectories in
->
+simulations. For instance, instantaneous temperature of an
->
+Hamiltonian system of $N$ particle can be measured by
->
+\[
->
+T = \sum\limits_{i = 1}^N {\frac{{m_i v_i^2 }}{{fk_B }}}
->
+\]
->
+where $m_i$ and $v_i$ are the mass and velocity of $i$th particle
->
+respectively, $f$ is the number of degrees of freedom, and $k_B$ is
->
+the boltzman constant.
-<
+\subsection{\label{introSec:forceEvaluation}Force Evaluation}
->
+A typical molecular dynamics run consists of three essential steps:
->
+\begin{enumerate}
->
+  \item Initialization
->
+    \begin{enumerate}
->
+    \item Preliminary preparation
->
+    \item Minimization
->
+    \item Heating
->
+    \item Equilibration
->
+    \end{enumerate}
->
+  \item Production
->
+  \item Analysis
->
+\end{enumerate}
->
+These three individual steps will be covered in the following
->
+sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
->
+initialization of a simulation. Sec.~\ref{introSec:production} will
->
+discusses issues in production run. Sec.~\ref{introSection:Analysis}
->
+provides the theoretical tools for trajectory analysis.
->
->
+\subsection{\label{introSec:initialSystemSettings}Initialization}
->
->
+\subsubsection{Preliminary preparation}
-<
+\subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}
->
+When selecting the starting structure of a molecule for molecular
->
+simulation, one may retrieve its Cartesian coordinates from public
->
+databases, such as RCSB Protein Data Bank \textit{etc}. Although
->
+thousands of crystal structures of molecules are discovered every
->
+year, many more remain unknown due to the difficulties of
->
+purification and crystallization. Even for the molecule with known
->
+structure, some important information is missing. For example, the
->
+missing hydrogen atom which acts as donor in hydrogen bonding must
->
+be added. Moreover, in order to include electrostatic interaction,
->
+one may need to specify the partial charges for individual atoms.
->
+Under some circumstances, we may even need to prepare the system in
->
+a special setup. For instance, when studying transport phenomenon in
->
+membrane system, we may prepare the lipids in bilayer structure
->
+instead of placing lipids randomly in solvent, since we are not
->
+interested in self-aggregation and it takes a long time to happen.
->
->
+\subsubsection{Minimization}
->
->
+It is quite possible that some of molecules in the system from
->
+preliminary preparation may be overlapped with each other. This
->
+close proximity leads to high potential energy which consequently
->
+jeopardizes any molecular dynamics simulations. To remove these
->
+steric overlaps, one typically performs energy minimization to find
->
+a more reasonable conformation. Several energy minimization methods
->
+have been developed to exploit the energy surface and to locate the
->
+local minimum. While converging slowly near the minimum, steepest
->
+descent method is extremely robust when systems are far from
->
+harmonic. Thus, it is often used to refine structure from
->
+crystallographic data. Relied on the gradient or hessian, advanced
->
+methods like conjugate gradient and Newton-Raphson converge rapidly
->
+to a local minimum, while become unstable if the energy surface is
->
+far from quadratic. Another factor must be taken into account, when
->
+choosing energy minimization method, is the size of the system.
->
+Steepest descent and conjugate gradient can deal with models of any
->
+size. Because of the limit of computation power to calculate hessian
->
+matrix and insufficient storage capacity to store them, most
->
+Newton-Raphson methods can not be used with very large models.
->
->
+\subsubsection{Heating}
-+
+Typically, Heating is performed by assigning random velocities
-+
+according to a Gaussian distribution for a temperature. Beginning at
-+
+a lower temperature and gradually increasing the temperature by
-+
+assigning greater random velocities, we end up with setting the
-+
+temperature of the system to a final temperature at which the
-+
+simulation will be conducted. In heating phase, we should also keep
-+
+the system from drifting or rotating as a whole. Equivalently, the
-+
+net linear momentum and angular momentum of the system should be
-+
+shifted to zero.
-+
-+
+\subsubsection{Equilibration}
-+
-+
+The purpose of equilibration is to allow the system to evolve
-+
+spontaneously for a period of time and reach equilibrium. The
-+
+procedure is continued until various statistical properties, such as
-+
+temperature, pressure, energy, volume and other structural
-+
+properties \textit{etc}, become independent of time. Strictly
-+
+speaking, minimization and heating are not necessary, provided the
-+
+equilibration process is long enough. However, these steps can serve
-+
+as a means to arrive at an equilibrated structure in an effective
-+
+way.
-+
-+
+\subsection{\label{introSection:production}Production}
-+
-+
+Production run is the most important steps of the simulation, in
-+
+which the equilibrated structure is used as a starting point and the
-+
+motions of the molecules are collected for later analysis. In order
-+
+to capture the macroscopic properties of the system, the molecular
-+
+dynamics simulation must be performed in correct and efficient way.
-+
-+
+The most expensive part of a molecular dynamics simulation is the
-+
+calculation of non-bonded forces, such as van der Waals force and
-+
+Coulombic forces \textit{etc}. For a system of $N$ particles, the
-+
+complexity of the algorithm for pair-wise interactions is $O(N^2 )$,
-+
+which making large simulations prohibitive in the absence of any
-+
+computation saving techniques.
-+
-+
+A natural approach to avoid system size issue is to represent the
-+
+bulk behavior by a finite number of the particles. However, this
-+
+approach will suffer from the surface effect. To offset this,
-+
+\textit{Periodic boundary condition} is developed to simulate bulk
-+
+properties with a relatively small number of particles. In this
-+
+method, the simulation box is replicated throughout space to form an
-+
+infinite lattice. During the simulation, when a particle moves in
-+
+the primary cell, its image in other cells move in exactly the same
-+
+direction with exactly the same orientation. Thus, as a particle
-+
+leaves the primary cell, one of its images will enter through the
-+
+opposite face.
-+
+%\begin{figure}
-+
+%\centering
-+
+%\includegraphics[width=\linewidth]{pbcFig.eps}
-+
+%\caption[An illustration of periodic boundary conditions]{A 2-D
-+
+%illustration of periodic boundary conditions. As one particle leaves
-+
+%the right of the simulation box, an image of it enters the left.}
-+
+%\label{introFig:pbc}
-+
+%\end{figure}
-+
-+
+%cutoff and minimum image convention
-+
+Another important technique to improve the efficiency of force
-+
+evaluation is to apply cutoff where particles farther than a
-+
+predetermined distance, are not included in the calculation
-+
+\cite{Frenkel1996}. The use of a cutoff radius will cause a
-+
+discontinuity in the potential energy curve
-+
+(Fig.~\ref{introFig:shiftPot}). Fortunately, one can shift the
-+
+potential to ensure the potential curve go smoothly to zero at the
-+
+cutoff radius. Cutoff strategy works pretty well for Lennard-Jones
-+
+interaction because of its short range nature. However, simply
-+
+truncating the electrostatic interaction with the use of cutoff has
-+
+been shown to lead to severe artifacts in simulations. Ewald
-+
+summation, in which the slowly conditionally convergent Coulomb
-+
+potential is transformed into direct and reciprocal sums with rapid
-+
+and absolute convergence, has proved to minimize the periodicity
-+
+artifacts in liquid simulations. Taking the advantages of the fast
-+
+Fourier transform (FFT) for calculating discrete Fourier transforms,
-+
+the particle mesh-based methods are accelerated from $O(N^{3/2})$ to
-+
+$O(N logN)$. An alternative approach is \emph{fast multipole
-+
+method}, which treats Coulombic interaction exactly at short range,
-+
+and approximate the potential at long range through multipolar
-+
+expansion. In spite of their wide acceptances at the molecular
-+
+simulation community, these two methods are hard to be implemented
-+
+correctly and efficiently. Instead, we use a damped and
-+
+charge-neutralized Coulomb potential method developed by Wolf and
-+
+his coworkers. The shifted Coulomb potential for particle $i$ and
-+
+particle $j$ at distance $r_{rj}$ is given by:
-+
+\begin{equation}
-+
+V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
-+
+r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
-+
+R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha
-+
+r_{ij})}{r_{ij}}\right\}. \label{introEquation:shiftedCoulomb}
-+
+\end{equation}
-+
+where $\alpha$ is the convergence parameter. Due to the lack of
-+
+inherent periodicity and rapid convergence,this method is extremely
-+
+efficient and easy to implement.
-+
+%\begin{figure}
-+
+%\centering
-+
+%\includegraphics[width=\linewidth]{pbcFig.eps}
-+
+%\caption[An illustration of shifted Coulomb potential]{An illustration of shifted Coulomb potential.}
-+
+%\label{introFigure:shiftedCoulomb}
-+
+%\end{figure}
-+
-+
+%multiple time step
-+
-+
+\subsection{\label{introSection:Analysis} Analysis}
-+
-+
+Recently, advanced visualization technique are widely applied to
-+
+monitor the motions of molecules. Although the dynamics of the
-+
+system can be described qualitatively from animation, quantitative
-+
+trajectory analysis are more appreciable. According to the
-+
+principles of Statistical Mechanics,
-+
+Sec.~\ref{introSection:statisticalMechanics}, one can compute
-+
+thermodynamics properties, analyze fluctuations of structural
-+
+parameters, and investigate time-dependent processes of the molecule
-+
+from the trajectories.
-+
-+
+\subsubsection{\label{introSection:thermodynamicsProperties}Thermodynamics Properties}
-+
-+
+Thermodynamics properties, which can be expressed in terms of some
-+
+function of the coordinates and momenta of all particles in the
-+
+system, can be directly computed from molecular dynamics. The usual
-+
+way to measure the pressure is based on virial theorem of Clausius
-+
+which states that the virial is equal to $-3Nk_BT$. For a system
-+
+with forces between particles, the total virial, $W$, contains the
-+
+contribution from external pressure and interaction between the
-+
+particles:
-+
+\[
-+
+W =  - 3PV + \left\langle {\sum\limits_{i < j} {r{}_{ij} \cdot
-+
+f_{ij} } } \right\rangle
-+
+\]
-+
+where $f_{ij}$ is the force between particle $i$ and $j$ at a
-+
+distance $r_{ij}$. Thus, the expression for the pressure is given
-+
+by:
-+
+\begin{equation}
-+
+P = \frac{{Nk_B T}}{V} - \frac{1}{{3V}}\left\langle {\sum\limits_{i
-+
+< j} {r{}_{ij} \cdot f_{ij} } } \right\rangle
-+
+\end{equation}
-+
-+
+\subsubsection{\label{introSection:structuralProperties}Structural Properties}
-+
-+
+Structural Properties of a simple fluid can be described by a set of
-+
+distribution functions. Among these functions,\emph{pair
-+
+distribution function}, also known as \emph{radial distribution
-+
+function}, is of most fundamental importance to liquid-state theory.
-+
+Pair distribution function can be gathered by Fourier transforming
-+
+raw data from a series of neutron diffraction experiments and
-+
+integrating over the surface factor \cite{Powles73}. The experiment
-+
+result can serve as a criterion to justify the correctness of the
-+
+theory. Moreover, various equilibrium thermodynamic and structural
-+
+properties can also be expressed in terms of radial distribution
-+
+function \cite{allen87:csl}.
-+
-+
+A pair distribution functions $g(r)$ gives the probability that a
-+
+particle $i$ will be located at a distance $r$ from a another
-+
+particle $j$ in the system
-+
+\[
-+
+g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
-+
+\ne i} {\delta (r - r_{ij} )} } } \right\rangle.
-+
+\]
-+
+Note that the delta function can be replaced by a histogram in
-+
+computer simulation. Figure
-+
+\ref{introFigure:pairDistributionFunction} shows a typical pair
-+
+distribution function for the liquid argon system. The occurrence of
-+
+several peaks in the plot of $g(r)$ suggests that it is more likely
-+
+to find particles at certain radial values than at others. This is a
-+
+result of the attractive interaction at such distances. Because of
-+
+the strong repulsive forces at short distance, the probability of
-+
+locating particles at distances less than about 2.5{\AA} from each
-+
+other is essentially zero.
-+
-+
+%\begin{figure}
-+
+%\centering
-+
+%\includegraphics[width=\linewidth]{pdf.eps}
-+
+%\caption[Pair distribution function for the liquid argon
-+
+%]{Pair distribution function for the liquid argon}
-+
+%\label{introFigure:pairDistributionFunction}
-+
+%\end{figure}
-+
-+
+\subsubsection{\label{introSection:timeDependentProperties}Time-dependent
-+
+Properties}
-+
-+
+Time-dependent properties are usually calculated using \emph{time
-+
+correlation function}, which correlates random variables $A$ and $B$
-+
+at two different time
-+
+\begin{equation}
-+
+C_{AB} (t) = \left\langle {A(t)B(0)} \right\rangle.
-+
+\label{introEquation:timeCorrelationFunction}
-+
+\end{equation}
-+
+If $A$ and $B$ refer to same variable, this kind of correlation
-+
+function is called \emph{auto correlation function}. One example of
-+
+auto correlation function is velocity auto-correlation function
-+
+which is directly related to transport properties of molecular
-+
+liquids:
-+
+\[
-+
+D = \frac{1}{3}\int\limits_0^\infty  {\left\langle {v(t) \cdot v(0)}
-+
+\right\rangle } dt
-+
+\]
-+
+where $D$ is diffusion constant. Unlike velocity autocorrelation
-+
+function which is averaging over time origins and over all the
-+
+atoms, dipole autocorrelation are calculated for the entire system.
-+
+The dipole autocorrelation function is given by:
-+
+\[
-+
+c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
-+
+\right\rangle
-+
+\]
-+
+Here $u_{tot}$ is the net dipole of the entire system and is given
-+
+by
-+
+\[
-+
+u_{tot} (t) = \sum\limits_i {u_i (t)}
-+
+\]
-+
+In principle, many time correlation functions can be related with
-+
+Fourier transforms of the infrared, Raman, and inelastic neutron
-+
+scattering spectra of molecular liquids. In practice, one can
-+
+extract the IR spectrum from the intensity of dipole fluctuation at
-+
+each frequency using the following relationship:
-+
+\[
-+
+\hat c_{dipole} (v) = \int_{ - \infty }^\infty  {c_{dipole} (t)e^{ -
-+
+i2\pi vt} dt}
-+
+\]
-+
+\section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
+Rigid bodies are frequently involved in the modeling of different
+The break through in geometric literature suggests that, in order to
+develop a long-term integration scheme, one should preserve the
+symplectic structure of the flow. Introducing conjugate momentum to
-<
+rotation matrix $A$ and re-formulating Hamiltonian's equation, a
->
+rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
+symplectic integrator, RSHAKE, was proposed to evolve the
+Hamiltonian system in a constraint manifold by iteratively
-<
+satisfying the orthogonality constraint $A_t A = 1$. An alternative
->
+satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
+method using quaternion representation was developed by Omelyan.
+However, both of these methods are iterative and inefficient. In
+this section, we will present a symplectic Lie-Poisson integrator
+& { - \sin \theta _1 } & {\cos \theta _1 }  \\
+\end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
+\]
-<
+To reduce the cost of computing expensive functions in e^{\Delta
-<
+tR_1 }, we can use Cayley transformation,
->
+To reduce the cost of computing expensive functions in $e^{\Delta
->
+tR_1 }$, we can use Cayley transformation,
+\[
+e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
+\]
-<
-<
+The flow maps for $T_2^r$ and $T_2^r$ can be found in the same
->
+The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
+manner.
+In order to construct a second-order symplectic method, we split the
+\varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
+_{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
+\]
-<
+Moreover, \varphi _{\Delta t/2,V} can be divided into two sub-flows
-<
+which corresponding to force and torque respectively,
->
+Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
->
+sub-flows which corresponding to force and torque respectively,
+\[
+\varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
+_{\Delta t/2,\tau }.
+\]
+Since the associated operators of $\varphi _{\Delta t/2,F} $ and
+$\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
-<
+order inside \varphi _{\Delta t/2,V} does not matter.
->
+order inside $\varphi _{\Delta t/2,V}$ does not matter.
+Furthermore, kinetic potential can be separated to translational
+kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
+mimics a simple heat bath with stochastic and dissipative forces,
+has been applied in a variety of studies. This section will review
+the theory of Langevin dynamics simulation. A brief derivation of
-<
+generalized Langevin Dynamics will be given first. Follow that, we
->
+generalized Langevin equation will be given first. Follow that, we
+will discuss the physical meaning of the terms appearing in the
+equation as well as the calculation of friction tensor from
+hydrodynamics theory.
-<
+\subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
->
+\subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
-+
+Harmonic bath model, in which an effective set of harmonic
-+
+oscillators are used to mimic the effect of a linearly responding
-+
+environment, has been widely used in quantum chemistry and
-+
+statistical mechanics. One of the successful applications of
-+
+Harmonic bath model is the derivation of Deriving Generalized
-+
+Langevin Dynamics. Lets consider a system, in which the degree of
-+
+freedom $x$ is assumed to couple to the bath linearly, giving a
-+
+Hamiltonian of the form
+\begin{equation}
+H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
-<
+\label{introEquation:bathGLE}
->
+\label{introEquation:bathGLE}.
+\end{equation}
-<
+where $H_B$ is harmonic bath Hamiltonian,
->
+Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
->
+with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
+\[
-<
+H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
-<
+}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
->
+H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
->
+}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
->
+\right\}}
+\]
-<
+and $\Delta U$ is bilinear system-bath coupling,
->
+where the index $\alpha$ runs over all the bath degrees of freedom,
->
+$\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
->
+the harmonic bath masses, and $\Delta U$ is bilinear system-bath
->
+coupling,
+\[
+\Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
+\]
-<
+Completing the square,
->
+where $g_\alpha$ are the coupling constants between the bath and the
->
+coordinate $x$. Introducing
+\[
-<
+H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
-<
+{\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
-<
+w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
-<
+w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
-<
+}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
-<
+\]
-<
+and putting it back into Eq.~\ref{introEquation:bathGLE},
->
+W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
->
+}}{{2m_\alpha  w_\alpha ^2 }}} x^2
->
+\] and combining the last two terms in Equation
->
+\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
->
+Hamiltonian as
+\[
+H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
+{\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
+w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
+w_\alpha ^2 }}x} \right)^2 } \right\}}
+\]
-–
+where
-–
+\[
-–
+W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
-–
+}}{{2m_\alpha  w_\alpha ^2 }}} x^2
-–
+\]
+Since the first two terms of the new Hamiltonian depend only on the
+system coordinates, we can get the equations of motion for
+Generalized Langevin Dynamics by Hamilton's equations
+\ref{introEquation:motionHamiltonianCoordinate,
+introEquation:motionHamiltonianMomentum},
-<
+\begin{align}
-<
+\dot p &=  - \frac{{\partial H}}{{\partial x}}
-<
+       &= m\ddot x
-<
+       &= - \frac{{\partial W(x)}}{{\partial x}} - \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)}
-<
+\label{introEquation:Lp5}
-<
+\end{align}
-<
+, and
-<
+\begin{align}
-<
+\dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
-<
+                &= m\ddot x_\alpha
-<
+                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
-<
+\end{align}
->
+\begin{equation}
->
+m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
->
+\sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
->
+\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)},
->
+\label{introEquation:coorMotionGLE}
->
+\end{equation}
->
+and
->
+\begin{equation}
->
+m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x_\alpha   -
->
+\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
->
+\label{introEquation:bathMotionGLE}
->
+\end{equation}
-<
+\subsection{\label{introSection:laplaceTransform}The Laplace Transform}
->
+In order to derive an equation for $x$, the dynamics of the bath
->
+variables $x_\alpha$ must be solved exactly first. As an integral
->
+transform which is particularly useful in solving linear ordinary
->
+differential equations, Laplace transform is the appropriate tool to
->
+solve this problem. The basic idea is to transform the difficult
->
+differential equations into simple algebra problems which can be
->
+solved easily. Then applying inverse Laplace transform, also known
->
+as the Bromwich integral, we can retrieve the solutions of the
->
+original problems.
-+
+Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
-+
+transform of f(t) is a new function defined as
+\[
-<
+L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
->
+L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
+\]
-+
+where  $p$ is real and  $L$ is called the Laplace Transform
-+
+Operator. Below are some important properties of Laplace transform
-+
+\begin{equation}
-+
+\begin{array}{c}
-+
+ L(x + y) = L(x) + L(y) \\
-+
+ L(ax) = aL(x) \\
-+
+ L(\dot x) = pL(x) - px(0) \\
-+
+ L(\ddot x) = p^2 L(x) - px(0) - \dot x(0) \\
-+
+ L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p) \\
-+
+ \end{array}
-+
+\end{equation}
-+
+Applying Laplace transform to the bath coordinates, we obtain
+\[
-<
+L(x + y) = L(x) + L(y)
->
+\begin{array}{c}
->
+ p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) \\
->
+ L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
->
+ \end{array}
+\]
-<
->
+By the same way, the system coordinates become
+\[
-<
+L(ax) = aL(x)
->
+\begin{array}{c}
->
+ mL(\ddot x) =  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
->
+  - \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x) - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0) - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}  \\
->
+ \end{array}
+\]
-+
+With the help of some relatively important inverse Laplace
-+
+transformations:
+\[
-<
+L(\dot x) = pL(x) - px(0)
->
+\begin{array}{c}
->
+ L(\cos at) = \frac{p}{{p^2  + a^2 }} \\
->
+ L(\sin at) = \frac{a}{{p^2  + a^2 }} \\
->
+ L(1) = \frac{1}{p} \\
->
+ \end{array}
+\]
-<
-<
+\[
-<
+L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
-<
+\]
-<
-<
+\[
-<
+L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
-<
+\]
-<
-<
+Some relatively important transformation,
-<
+\[
-<
+L(\cos at) = \frac{p}{{p^2  + a^2 }}
-<
+\]
-<
-<
+\[
-<
+L(\sin at) = \frac{a}{{p^2  + a^2 }}
-<
+\]
-<
-<
+\[
-<
+L(1) = \frac{1}{p}
-<
+\]
-<
-<
+First, the bath coordinates,
-<
+\[
-<
+p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
-<
+_\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
-<
+}}L(x)
-<
+\]
-<
+\[
-<
+L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
-<
+px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
-<
+\]
-<
+Then, the system coordinates,
-<
+\begin{align}
-<
+mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
-<
+\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
-<
+}}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
-<
+(0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
-<
+}}\omega _\alpha ^2 L(x)} \right\}}
-<
+%
-<
+ &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
-<
+ \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
-<
+ - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
-<
+ - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
-<
+\end{align}
-<
+Then, the inverse transform,
-<
->
+, we obtain
+\begin{align}
+m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
+\sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
+(\omega _\alpha  t)} \right\}}
+\end{align}
-+
+Introducing a \emph{dynamic friction kernel}
+\begin{equation}
-+
+\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
-+
+}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
-+
+\label{introEquation:dynamicFrictionKernelDefinition}
-+
+\end{equation}
-+
+and \emph{a random force}
-+
+\begin{equation}
-+
+R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
-+
+- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
-+
+\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
-+
+(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t),
-+
+\label{introEquation:randomForceDefinition}
-+
+\end{equation}
-+
+the equation of motion can be rewritten as
-+
+\begin{equation}
+m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
+(t)\dot x(t - \tau )d\tau }  + R(t)
+\label{introEuqation:GeneralizedLangevinDynamics}
+\end{equation}
-<
+%where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
-<
+%$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
->
+which is known as the \emph{generalized Langevin equation}.
->
->
+\subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
->
->
+One may notice that $R(t)$ depends only on initial conditions, which
->
+implies it is completely deterministic within the context of a
->
+harmonic bath. However, it is easy to verify that $R(t)$ is totally
->
+uncorrelated to $x$ and $\dot x$,
+\[
-<
+\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
-<
+}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
->
+\begin{array}{l}
->
+ \left\langle {x(t)R(t)} \right\rangle  = 0, \\
->
+ \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
->
+ \end{array}
+\]
-<
+For an infinite harmonic bath, we can use the spectral density and
-<
+an integral over frequencies.
->
+This property is what we expect from a truly random process. As long
->
+as the model, which is gaussian distribution in general, chosen for
->
+$R(t)$ is a truly random process, the stochastic nature of the GLE
->
+still remains.
-+
+%dynamic friction kernel
-+
+The convolution integral
+\[
-<
+R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
-<
+- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
-<
+\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
-<
+(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
->
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }
+\]
-<
+The random forces depend only on initial conditions.
->
+depends on the entire history of the evolution of $x$, which implies
->
+that the bath retains memory of previous motions. In other words,
->
+the bath requires a finite time to respond to change in the motion
->
+of the system. For a sluggish bath which responds slowly to changes
->
+in the system coordinate, we may regard $\xi(t)$ as a constant
->
+$\xi(t) = \Xi_0$. Hence, the convolution integral becomes
->
+\[
->
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
->
+\]
->
+and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
->
+\[
->
+m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
->
+\frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
->
+\]
->
+which can be used to describe dynamic caging effect. The other
->
+extreme is the bath that responds infinitely quickly to motions in
->
+the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
->
+time:
->
+\[
->
+\xi (t) = 2\xi _0 \delta (t)
->
+\]
->
+Hence, the convolution integral becomes
->
+\[
->
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
->
+{\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
->
+\]
->
+and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
->
+\begin{equation}
->
+m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
->
+x(t) + R(t) \label{introEquation:LangevinEquation}
->
+\end{equation}
->
+which is known as the Langevin equation. The static friction
->
+coefficient $\xi _0$ can either be calculated from spectral density
->
+or be determined by Stokes' law for regular shaped particles.A
->
+briefly review on calculating friction tensor for arbitrary shaped
->
+particles is given in Sec.~\ref{introSection:frictionTensor}.
+\subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
-<
+So we can define a new set of coordinates,
->
->
+Defining a new set of coordinates,
+\[
+q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
+^2 }}x(0)
-<
+\]
-<
+This makes
->
+\],
->
+we can rewrite $R(T)$ as
+\[
-<
+R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
->
+R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
+\]
+And since the $q$ coordinates are harmonic oscillators,
+\[
-<
+\begin{array}{l}
->
+\begin{array}{c}
->
+ \left\langle {q_\alpha ^2 } \right\rangle  = \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
+ \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
+ \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
-+
+ \left\langle {R(t)R(0)} \right\rangle  = \sum\limits_\alpha  {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle } }  \\
-+
+  = \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
-+
+  = kT\xi (t) \\
+ \end{array}
+\]
-<
-<
+\begin{align}
-<
+\left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
-<
+{\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
-<
+(t)q_\beta  (0)} \right\rangle } }
-<
+%
-<
+&= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
-<
+\right\rangle \cos (\omega _\alpha  t)}
-<
+%
-<
+&= kT\xi (t)
-<
+\end{align}
-<
->
+Thus, we recover the \emph{second fluctuation dissipation theorem}
+\begin{equation}
+\xi (t) = \left\langle {R(t)R(0)} \right\rangle
-<
+\label{introEquation:secondFluctuationDissipation}
->
+\label{introEquation:secondFluctuationDissipation}.
+\end{equation}
-+
+In effect, it acts as a constraint on the possible ways in which one
-+
+can model the random force and friction kernel.
+\subsection{\label{introSection:frictionTensor} Friction Tensor}
+Theoretically, the friction kernel can be determined using velocity
+B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
+)T_{ij}
+\]
-<
+where \delta _{ij} is Kronecker delta function. Inverting matrix
->
+where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
+$B$, we obtain
+\[
+Form 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 do depend on the
->
+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
+obtain the resistance tensor at $P$ by
+\]
+where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
+joining center of resistance $R$ and origin $O$.
-–
-–
+%\section{\label{introSection:correlationFunctions}Correlation Functions}

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing trunk/tengDissertation/Introduction.tex (file contents): Revision 2718 by tim, Tue Apr 18 04:11:56 2006 UTC vs. Revision 2725 by tim, Fri Apr 21 05:45:14 2006 UTC

Diff Legend

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2718 by tim, Tue Apr 18 04:11:56 2006 UTC vs.
Revision 2725 by tim, Fri Apr 21 05:45:14 2006 UTC