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

Comparing trunk/xDissertation/Introduction.tex (file contents):
Revision 3336 by xsun, Wed Jan 30 16:01:02 2008 UTC vs.
Revision 3372 by xsun, Wed Mar 19 21:04:47 2008 UTC

+\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
-+
-+
+\section{Background on the Problem\label{In:sec:pro}}
-+
+Phospholipid molecules are the primary topic of this dissertation
-+
+because of their critical role as the foundation of biological
-+
+membranes. Lipids, when dispersed in water, self assemble into a
-+
+mumber of topologically distinct bilayer structures. The phase
-+
+behavior of lipid bilayers has been explored
-+
+experimentally~\cite{Cevc87}, however, a complete understanding of the
-+
+mechanism and driving forces behind the various phases has not been
-+
+achieved.
-+
-+
+\subsection{Ripple Phase\label{In:ssec:ripple}}
-+
+The $P_{\beta'}$ {\it ripple phase} of lipid bilayers, named from the
-+
+periodic buckling of the membrane, is an intermediate phase which is
-+
+developed either from heating the gel phase $L_{\beta'}$ or cooling
-+
+the fluid phase $L_\alpha$. A Sketch is shown in
-+
+figure~\ref{Infig:phaseDiagram}.~\cite{Cevc87}
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=\linewidth]{./figures/inPhaseDiagram.pdf}
-+
+\caption{A phase diagram of lipid bilayer. With increasing the
-+
+temperature, the bilayer can go through a gel ($L_{\beta'}$), ripple
-+
+($P_{\beta'}$) to fluid ($L_\alpha$) phase transition.}
-+
+\label{Infig:phaseDiagram}
-+
+\end{figure}
-+
+Most structural information of the ripple phase has been obtained by
-+
+the X-ray diffraction~\cite{Sun96,Katsaras00} and freeze-fracture
-+
+electron microscopy (FFEM).~\cite{Copeland80,Meyer96} The X-ray
-+
+diffraction work by Katsaras {\it et al.} showed that a rich phase
-+
+diagram exhibiting both {\it asymmetric} and {\it symmetric} ripples
-+
+is possible for lecithin bilayers.\cite{Katsaras00} Recently,
-+
+Kaasgaard {\it et al.} used atomic force microscopy (AFM) to observe
-+
+ripple phase morphology in bilayers supported on
-+
+mica.~\cite{Kaasgaard03}
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=\linewidth]{./figures/inRipple.pdf}
-+
+\caption{The experimental observed ripple phase. The top image is
-+
+obtained by X-ray diffraction~\cite{Sun96}, and the bottom one is
-+
+observed by AFM.~\cite{Kaasgaard03}}
-+
+\label{Infig:ripple}
-+
+\end{figure}
-+
+Figure~\ref{Infig:ripple} shows the ripple phase oberved by X-ray
-+
+diffraction and AFM. The experimental results provide strong support
-+
+for a 2-dimensional triangular packing lattice of the lipid molecules
-+
+within the ripple phase.  This is a notable change from the observed
-+
+lipid packing within the gel phase,~\cite{Cevc87} although Tenchov
-+
+{\it et al.} have recently observed near-hexagonal packing in some
-+
+phosphatidylcholine (PC) gel phases.~\cite{Tenchov2001} However, the
-+
+physical mechanism for the formation of the ripple phase has never
-+
+been explained and the microscopic structure of the ripple phase has
-+
+never been elucidated by experiments. Computational simulation is a
-+
+perfect tool to study the microscopic properties for a
-+
+system. However, the large length scale the ripple structure and the
-+
+long time scale of the formation of the ripples are crucial obstacles
-+
+to performing the actual work. The principal ideas explored in this
-+
+dissertation are attempts to break the computational task up by
-+
+\begin{itemize}
-+
+\item Simplifying the lipid model.
-+
+\item Improving algorithm for integrating the equations of motion.
-+
+\end{itemize}
-+
+In chapters~\ref{chap:mc} and~\ref{chap:md}, we use a simple point
-+
+dipole spin model and a coarse-grained molecualr scale model to
-+
+perform the Monte Carlo and Molecular Dynamics simulations
-+
+respectively, and in chapter~\ref{chap:ld}, we develop a Langevin
-+
+Dynamics algorithm which excludes the explicit solvent to improve the
-+
+efficiency of the simulations.
-+
-+
+\subsection{Lattice Model\label{In:ssec:model}}
-+
+The gel-like characteristic (relatively immobile molecules) exhibited
-+
+by the ripple phase makes it feasible to apply a lattice model to
-+
+study the system. The popular $2$ dimensional lattice models, {\it
-+
+i.e.}, the Ising, $X-Y$, and Heisenberg models, show {\it frustration}
-+
+on triangular lattices. The Hamiltonians of these systems are given by
-+
+\begin{equation}
-+
+H =
-+
+  \begin{cases}
-+
+    -J \sum_n \sum_{n'} s_n s_n' & \text{Ising}, \\
-+
+    -J \sum_n \sum_{n'} \vec s_n \cdot \vec s_{n'} & \text{$X-Y$ and
-+
+Heisenberg},
-+
+  \end{cases}
-+
+\end{equation}
-+
+where $J$ has non zero value only when spins $s_n$ ($\vec s_n$) and
-+
+$s_{n'}$ ($\vec s_{n'}$) are nearest neighbors. When $J > 0$, spins
-+
+prefer an aligned structure, and if $J < 0$, spins prefer an
-+
+anti-aligned structure.
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=3in]{./figures/inFrustration.pdf}
-+
+\caption{Frustration on triangular lattice, the spins and dipoles are
-+
+represented by arrows. The multiple local minima of energy states
-+
+induce the frustration for spins and dipoles picking the directions.}
-+
+\label{Infig:frustration}
-+
+\end{figure}
-+
+The spins in figure~\ref{Infig:frustration} shows an illustration of
-+
+the frustration for $J < 0$ on a triangular lattice. There are
-+
+multiple local minima energy states which are independent of the
-+
+direction of the spin on top of the triangle, therefore infinite
-+
+possibilities for the packing of spins which induces what is known as
-+
+``complete regular frustration'' which leads to disordered low
-+
+temperature phases. The similarity goes to the dipoles on a hexagonal
-+
+lattice, which are shown by the dipoles in
-+
+figure~\ref{Infig:frustration}. In this circumstance, the dipoles want
-+
+to be aligned, however, due to the long wave fluctuation, at low
-+
+temperature, the aligned state becomes unstable, vortex is formed and
-+
+results in multiple local minima of energy states. The dipole on the
-+
+center of the hexagonal lattice is frustrated.
-+
-+
+The lack of translational degrees of freedom in lattice models
-+
+prevents their utilization in models for surface buckling. In
-+
+chapter~\ref{chap:mc}, a modified lattice model is introduced to
-+
+tackle this specific situation.
-+
-+
+\section{Overview of Classical Statistical Mechanics\label{In:sec:SM}}
-+
+Statistical mechanics provides a way to calculate the macroscopic
-+
+properties of a system from the molecular interactions used in
-+
+computational simulations. This section serves as a brief introduction
-+
+to key concepts of classical statistical mechanics that we used in
-+
+this dissertation. Tolman gives an excellent introduction to the
-+
+principles of statistical mechanics.~\cite{Tolman1979} A large part of
-+
+section~\ref{In:sec:SM} will follow Tolman's notation.
-+
-+
+\subsection{Ensembles\label{In:ssec:ensemble}}
-+
+In classical mechanics, the state of the system is completely
-+
+described by the positions and momenta of all particles. If we have an
-+
+$N$ particle system, there are $6N$ coordinates ($3N$ positions $(q_1,
-+
+q_2, \ldots, q_{3N})$ and $3N$ momenta $(p_1, p_2, \ldots, p_{3N})$)
-+
+to define the instantaneous state of the system. Each single set of
-+
+the $6N$ coordinates can be considered as a unique point in a $6N$
-+
+dimensional space where each perpendicular axis is one of
-+
+$q_{i\alpha}$ or $p_{i\alpha}$ ($i$ is the particle and $\alpha$ is
-+
+the spatial axis). This $6N$ dimensional space is known as phase
-+
+space. The instantaneous state of the system is a single point in
-+
+phase space. A trajectory is a connected path of points in phase
-+
+space. An ensemble is a collection of systems described by the same
-+
+macroscopic observables but which have microscopic state
-+
+distributions. In phase space an ensemble is a collection of a set of
-+
+representive points. A density distribution $\rho(q^N, p^N)$ of the
-+
+representive points in phase space describes the condition of an
-+
+ensemble of identical systems. Since this density may change with
-+
+time, it is also a function of time. $\rho(q^N, p^N, t)$ describes the
-+
+ensemble at a time $t$, and $\rho(q^N, p^N, t')$ describes the
-+
+ensemble at a later time $t'$. For convenience, we will use $\rho$
-+
+instead of $\rho(q^N, p^N, t)$ in the following disccusion. If we
-+
+normalize $\rho$ to unity,
-+
+\begin{equation}
-+
+= \int d \vec q~^N \int d \vec p~^N \rho,
-+
+\label{Ineq:normalized}
-+
+\end{equation}
-+
+then the value of $\rho$ gives the probability of finding the system
-+
+in a unit volume in phase space.
-+
-+
+Liouville's theorem describes the change in density $\rho$ with
-+
+time. The number of representative points at a given volume in phase
-+
+space at any instant can be written as:
-+
+\begin{equation}
-+
+\label{Ineq:deltaN}
-+
+\delta N = \rho~\delta q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N.
-+
+\end{equation}
-+
+To calculate the change in the number of representative points in this
-+
+volume, let us consider a simple condition: the change in the number
-+
+of representative points along the $q_1$ axis. The rate of the number
-+
+of the representative points entering the volume at $q_1$ per unit time
-+
+is:
-+
+\begin{equation}
-+
+\label{Ineq:deltaNatq1}
-+
+\rho~\dot q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N,
-+
+\end{equation}
-+
+and the rate of the number of representative points leaving the volume
-+
+at another position $q_1 + \delta q_1$ is:
-+
+\begin{equation}
-+
+\label{Ineq:deltaNatq1plusdeltaq1}
-+
+\left( \rho + \frac{\partial \rho}{\partial q_1} \delta q_1 \right)\left(\dot q_1 +
-+
+\frac{\partial \dot q_1}{\partial q_1} \delta q_1 \right)\delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N.
-+
+\end{equation}
-+
+Here the higher order differentials are neglected. So the change in
-+
+the number of representative points is the difference between
-+
+eq.~\ref{Ineq:deltaNatq1} and eq.~\ref{Ineq:deltaNatq1plusdeltaq1},
-+
+which gives us:
-+
+\begin{equation}
-+
+\label{Ineq:deltaNatq1axis}
-+
+-\left(\rho \frac{\partial {\dot q_1}}{\partial q_1} + \frac{\partial {\rho}}{\partial q_1} \dot q_1 \right)\delta q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N,
-+
+\end{equation}
-+
+where, higher order differetials are neglected. If we sum over all the
-+
+axes in the phase space, we can get the change in the number of
-+
+representative points in a given volume with time:
-+
+\begin{equation}
-+
+\label{Ineq:deltaNatGivenVolume}
-+
+\frac{d(\delta N)}{dt} = -\sum_{i=1}^N \left[\rho \left(\frac{\partial
-+
+{\dot q_i}}{\partial q_i} + \frac{\partial
-+
+{\dot p_i}}{\partial p_i}\right) + \left( \frac{\partial {\rho}}{\partial
-+
+q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i\right)\right]\delta q_1 \delta q_2 \ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N.
-+
+\end{equation}
-+
+From Hamilton's equations of motion,
-+
+\begin{equation}
-+
+\frac{\partial {\dot q_i}}{\partial q_i} = - \frac{\partial
-+
+{\dot p_i}}{\partial p_i},
-+
+\label{Ineq:canonicalFormOfEquationOfMotion}
-+
+\end{equation}
-+
+this cancels out the first term on the right side of
-+
+eq.~\ref{Ineq:deltaNatGivenVolume}. If both sides of
-+
+eq.~\ref{Ineq:deltaNatGivenVolume} are divided by $\delta q_1 \delta q_2
-+
+\ldots \delta q_N \delta p_1 \delta p_2 \ldots \delta p_N$, then we
-+
+arrive at Liouville's theorem:
-+
+\begin{equation}
-+
+\left( \frac{\partial \rho}{\partial t} \right)_{q, p} = -\sum_{i} \left(
-+
+\frac{\partial {\rho}}{\partial
-+
+q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i \right).
-+
+\label{Ineq:simpleFormofLiouville}
-+
+\end{equation}
-+
+This is the basis of statistical mechanics. If we move the right
-+
+side of equation~\ref{Ineq:simpleFormofLiouville} to the left, we
-+
+will obtain
-+
+\begin{equation}
-+
+\left( \frac{\partial \rho}{\partial t} \right)_{q, p} + \sum_{i} \left(
-+
+\frac{\partial {\rho}}{\partial
-+
+q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i \right)
-+
+= 0.
-+
+\label{Ineq:anotherFormofLiouville}
-+
+\end{equation}
-+
+It is easy to note that the left side of
-+
+equation~\ref{Ineq:anotherFormofLiouville} is the total derivative of
-+
+$\rho$ with respect to $t$, which means
-+
+\begin{equation}
-+
+\frac{d \rho}{dt} = 0,
-+
+\label{Ineq:conservationofRho}
-+
+\end{equation}
-+
+and the rate of density change is zero in the neighborhood of any
-+
+selected moving representative points in the phase space.
-+
-+
+The condition of the ensemble is determined by the density
-+
+distribution. If we consider the density distribution as only a
-+
+function of $q$ and $p$, which means the rate of change of the phase
-+
+space density in the neighborhood of all representative points in the
-+
+phase space is zero,
-+
+\begin{equation}
-+
+\left( \frac{\partial \rho}{\partial t} \right)_{q, p} = 0.
-+
+\label{Ineq:statEquilibrium}
-+
+\end{equation}
-+
+We may conclude the ensemble is in {\it statistical equilibrium}. An
-+
+ensemble in statistical equilibrium means the system is also in
-+
+macroscopic equilibrium. If $\left( \frac{\partial \rho}{\partial t}
-+
+\right)_{q, p} = 0$, then
-+
+\begin{equation}
-+
+ \sum_{i} \left(
-+
+\frac{\partial {\rho}}{\partial
-+
+q_i} \dot q_i + \frac{\partial {\rho}}{\partial p_i} \dot p_i \right)
-+
+= 0.
-+
+\label{Ineq:constantofMotion}
-+
+\end{equation}
-+
+If $\rho$ is a function only of some constant of the motion, $\rho$ is
-+
+independent of time. For a conservative system, the energy of the
-+
+system is one of the constants of the motion. There are many
-+
+thermodynamically relevant ensembles: when the density distribution is
-+
+constant everywhere in the phase space,
-+
+\begin{equation}
-+
+\rho = \mathrm{const.}
-+
+\label{Ineq:uniformEnsemble}
-+
+\end{equation}
-+
+the ensemble is called {\it uniform ensemble}.
-+
-+
+\subsubsection{The Microcanonical Ensemble\label{In:sssec:microcanonical}}
-+
+Another useful ensemble is the {\it microcanonical ensemble}, for
-+
+which:
-+
+\begin{equation}
-+
+\rho = \delta \left( H(q^N, p^N) - E \right) \frac{1}{\Sigma (N, V, E)}
-+
+\label{Ineq:microcanonicalEnsemble}
-+
+\end{equation}
-+
+where $\Sigma(N, V, E)$ is a normalization constant parameterized by
-+
+$N$, the total number of particles, $V$, the total physical volume and
-+
+$E$, the total energy. The physical meaning of $\Sigma(N, V, E)$ is
-+
+the phase space volume accessible to a microcanonical system with
-+
+energy $E$ evolving under Hamilton's equations. $H(q^N, p^N)$ is the
-+
+Hamiltonian of the system. The Gibbs entropy is defined as
-+
+\begin{equation}
-+
+S = - k_B \int d \vec q~^N \int d \vec p~^N \rho \ln [C^N \rho],
-+
+\label{Ineq:gibbsEntropy}
-+
+\end{equation}
-+
+where $k_B$ is the Boltzmann constant and $C^N$ is a number which
-+
+makes the argument of $\ln$ dimensionless. In this case, $C^N$ is the
-+
+total phase space volume of one state. The entropy of a microcanonical
-+
+ensemble is given by
-+
+\begin{equation}
-+
+S = k_B \ln \left(\frac{\Sigma(N, V, E)}{C^N}\right).
-+
+\label{Ineq:entropy}
-+
+\end{equation}
-+
-+
+\subsubsection{The Canonical Ensemble\label{In:sssec:canonical}}
-+
+If the density distribution $\rho$ is given by
-+
+\begin{equation}
-+
+\rho = \frac{1}{Z_N}e^{-H(q^N, p^N) / k_B T},
-+
+\label{Ineq:canonicalEnsemble}
-+
+\end{equation}
-+
+the ensemble is known as the {\it canonical ensemble}. Here,
-+
+\begin{equation}
-+
+Z_N = \int d \vec q~^N \int_\Gamma d \vec p~^N  e^{-H(q^N, p^N) / k_B T},
-+
+\label{Ineq:partitionFunction}
-+
+\end{equation}
-+
+which is also known as the canonical{\it partition function}. $\Gamma$
-+
+indicates that the integral is over all phase space. In the canonical
-+
+ensemble, $N$, the total number of particles, $V$, total volume, and
-+
+$T$, the temperature, are constants. The systems with the lowest
-+
+energies hold the largest population. According to maximum principle,
-+
+thermodynamics maximizes the entropy $S$, implying that
-+
+\begin{equation}
-+
+\begin{array}{ccc}
-+
+\delta S = 0 & \mathrm{and} & \delta^2 S < 0.
-+
+\end{array}
-+
+\label{Ineq:maximumPrinciple}
-+
+\end{equation}
-+
+From Eq.~\ref{Ineq:maximumPrinciple} and two constrains on the
-+
+canonical ensemble, {\it i.e.}, total probability and average energy
-+
+must be conserved, the partition function can be shown to be
-+
+equivalent to
-+
+\begin{equation}
-+
+Z_N = e^{-A/k_B T},
-+
+\label{Ineq:partitionFunctionWithFreeEnergy}
-+
+\end{equation}
-+
+where $A$ is the Helmholtz free energy. The significance of
-+
+Eq.~\ref{Ineq:entropy} and~\ref{Ineq:partitionFunctionWithFreeEnergy} is
-+
+that they serve as a connection between macroscopic properties of the
-+
+system and the distribution of microscopic states.
-+
-+
+There is an implicit assumption that our arguments are based on so
-+
+far. A representative point in the phase space is equally likely to be
-+
+found in any energetically allowed region. In other words, all
-+
+energetically accessible states are represented equally, the
-+
+probabilities to find the system in any of the accessible states is
-+
+equal. This is called the principle of equal a {\it priori}
-+
+probabilities.
-+
-+
+\subsection{Statistical Averages\label{In:ssec:average}}
-+
+Given a density distribution $\rho$ in phase space, the average
-+
+of any quantity ($F(q^N, p^N$)) which depends on the coordinates
-+
+($q^N$) and the momenta ($p^N$) for all the systems in the ensemble
-+
+can be calculated based on the definition shown by
-+
+Eq.~\ref{Ineq:statAverage1}
-+
+\begin{equation}
-+
+\langle F(q^N, p^N) \rangle = \frac{\int d \vec q~^N \int d \vec p~^N
-+
+F(q^N, p^N) \rho}{\int d \vec q~^N \int d \vec p~^N \rho}.
-+
+\label{Ineq:statAverage1}
-+
+\end{equation}
-+
+Since the density distribution $\rho$ is normalized to unity, the mean
-+
+value of $F(q^N, p^N)$ is simplified to
-+
+\begin{equation}
-+
+\langle F(q^N, p^N) \rangle = \int d \vec q~^N \int d \vec p~^N F(q^N,
-+
+p^N) \rho,
-+
+\label{Ineq:statAverage2}
-+
+\end{equation}
-+
+called the {\it ensemble average}. However, the quantity is often
-+
+averaged for a finite time in real experiments,
-+
+\begin{equation}
-+
+\langle F(q^N, p^N) \rangle_t = \lim_{T \rightarrow \infty}
-+
+\frac{1}{T} \int_{t_0}^{t_0+T} F[q^N(t), p^N(t)] dt.
-+
+\label{Ineq:timeAverage1}
-+
+\end{equation}
-+
+Usually this time average is independent of $t_0$ in statistical
-+
+mechanics, so Eq.~\ref{Ineq:timeAverage1} becomes
-+
+\begin{equation}
-+
+\langle F(q^N, p^N) \rangle_t = \lim_{T \rightarrow \infty}
-+
+\frac{1}{T} \int_{0}^{T} F[q^N(t), p^N(t)] dt
-+
+\label{Ineq:timeAverage2}
-+
+\end{equation}
-+
+for an infinite time interval.
-+
-+
+\subsubsection{Ergodicity\label{In:sssec:ergodicity}}
-+
+The {\it ergodic hypothesis}, an important hypothesis governing modern
-+
+computer simulations states that the system will eventually pass
-+
+arbitrarily close to any point that is energetically accessible in
-+
+phase space. Mathematically, this leads to
-+
+\begin{equation}
-+
+\langle F(q^N, p^N) \rangle = \langle F(q^N, p^N) \rangle_t.
-+
+\label{Ineq:ergodicity}
-+
+\end{equation}
-+
+Eq.~\ref{Ineq:ergodicity} validates Molecular Dynamics as a form of
-+
+averaging for sufficiently ergodic systems. Also Monte Carlo may be
-+
+used to obtain ensemble averages of a quantity which are related to
-+
+time averages measured in experiments.
-+
-+
+\subsection{Correlation Functions\label{In:ssec:corr}}
-+
+Thermodynamic properties can be computed by equilibrium statistical
-+
+mechanics. On the other hand, {\it Time correlation functions} are a
-+
+powerful tool to understand the evolution of a dynamical
-+
+systems. Imagine that property $A(q^N, p^N, t)$ as a function of
-+
+coordinates $q^N$ and momenta $p^N$ has an intial value at $t_0$, and
-+
+at a later time $t_0 + \tau$ this value has changed. If $\tau$ is very
-+
+small, the change of the value is minor, and the later value of
-+
+$A(q^N, p^N, t_0 + \tau)$ is correlated to its initial value. However,
-+
+when $\tau$ is large, this correlation is lost. A time correlation
-+
+function measures this relationship and is defined
-+
+by~\cite{Berne90}
-+
+\begin{equation}
-+
+C_{AA}(\tau) = \langle A(0)A(\tau) \rangle = \lim_{T \rightarrow
-+
+\infty}
-+
+\frac{1}{T} \int_{0}^{T} dt A(t) A(t + \tau).
-+
+\label{Ineq:autocorrelationFunction}
-+
+\end{equation}
-+
+Eq.~\ref{Ineq:autocorrelationFunction} is the correlation function of
-+
+a single variable, called an {\it autocorrelation function}. The
-+
+definition of the correlation function for two different variables is
-+
+similar to that of autocorrelation function, which is
-+
+\begin{equation}
-+
+C_{AB}(\tau) = \langle A(0)B(\tau) \rangle = \lim_{T \rightarrow \infty}
-+
+\frac{1}{T} \int_{0}^{T} dt A(t) B(t + \tau),
-+
+\label{Ineq:crosscorrelationFunction}
-+
+\end{equation}
-+
+and called {\it cross correlation function}.
-+
-+
+We know from the ergodic hypothesis that there is a relationship
-+
+between time average and ensemble average. We can put the correlation
-+
+function in a classical mechanics form,
-+
+\begin{equation}
-+
+C_{AA}(\tau) = \int d \vec q~^N \int d \vec p~^N A[(q^N(t), p^N(t)]
-+
+A[q^N(t+\tau), q^N(t+\tau)] \rho(q, p)
-+
+\label{Ineq:autocorrelationFunctionCM}
-+
+\end{equation}
-+
+and
-+
+\begin{equation}
-+
+C_{AB}(\tau) = \int d \vec q~^N \int d \vec p~^N A[(q^N(t), p^N(t)]
-+
+B[q^N(t+\tau), q^N(t+\tau)] \rho(q, p)
-+
+\label{Ineq:crosscorrelationFunctionCM}
-+
+\end{equation}
-+
+as autocorrelation function and cross correlation function
-+
+respectively. $\rho(q, p)$ is the density distribution at equillibrium
-+
+in phase space. In many cases, the correlation function decay is a
-+
+single exponential
-+
+\begin{equation}
-+
+C(t) \sim e^{-t / \tau_r},
-+
+\label{Ineq:relaxation}
-+
+\end{equation}
-+
+where $\tau_r$ is known as relaxation time which discribes the rate of
-+
+the decay.
-+
-+
+\section{Methodolody\label{In:sec:method}}
-+
+The simulations performed in this dissertation are branched into two
-+
+main catalog, Monte Carlo and Molecular Dynamics. There are two main
-+
+difference between Monte Carlo and Molecular Dynamics simulations. One
-+
+is that the Monte Carlo simulation is time independent, and Molecular
-+
+Dynamics simulation is time involved. Another dissimilar is that the
-+
+Monte Carlo is a stochastic process, the configuration of the system
-+
+is not determinated by its past, however, using Moleuclar Dynamics,
-+
+the system is propagated by Newton's equation of motion, the
-+
+trajectory of the system evolved in the phase space is determined. A
-+
+brief introduction of the two algorithms are given in
-+
+section~\ref{In:ssec:mc} and~\ref{In:ssec:md}. An extension of the
-+
+Molecular Dynamics, Langevin Dynamics, is introduced by
-+
+section~\ref{In:ssec:ld}.
-+
-+
+\subsection{Monte Carlo\label{In:ssec:mc}}
-+
+Monte Carlo algorithm was first introduced by Metropolis {\it et
-+
+al.}.~\cite{Metropolis53} Basic Monte Carlo algorithm is usually
-+
+applied to the canonical ensemble, a Boltzmann-weighted ensemble, in
-+
+which the $N$, the total number of particles, $V$, total volume, $T$,
-+
+temperature are constants. The average energy is given by substituding
-+
+Eq.~\ref{Ineq:canonicalEnsemble} into Eq.~\ref{Ineq:statAverage2},
-+
+\begin{equation}
-+
+\langle E \rangle = \frac{1}{Z_N} \int d \vec q~^N \int d \vec p~^N E e^{-H(q^N, p^N) / k_B T}.
-+
+\label{Ineq:energyofCanonicalEnsemble}
-+
+\end{equation}
-+
+So are the other properties of the system. The Hamiltonian is the
-+
+summation of Kinetic energy $K(p^N)$ as a function of momenta and
-+
+Potential energy $U(q^N)$ as a function of positions,
-+
+\begin{equation}
-+
+H(q^N, p^N) = K(p^N) + U(q^N).
-+
+\label{Ineq:hamiltonian}
-+
+\end{equation}
-+
+If the property $A$ is only a function of position ($ A = A(q^N)$),
-+
+the mean value of $A$ is reduced to
-+
+\begin{equation}
-+
+\langle A \rangle = \frac{\int d \vec q~^N \int d \vec p~^N A e^{-U(q^N) / k_B T}}{\int d \vec q~^N \int d \vec p~^N e^{-U(q^N) / k_B T}},
-+
+\label{Ineq:configurationIntegral}
-+
+\end{equation}
-+
+The kinetic energy $K(p^N)$ is factored out in
-+
+Eq.~\ref{Ineq:configurationIntegral}. $\langle A
-+
+\rangle$ is a configuration integral now, and the
-+
+Eq.~\ref{Ineq:configurationIntegral} is equivalent to
-+
+\begin{equation}
-+
+\langle A \rangle = \int d \vec q~^N A \rho(q^N).
-+
+\label{Ineq:configurationAve}
-+
+\end{equation}
-+
-+
+In a Monte Carlo simulation of canonical ensemble, the probability of
-+
+the system being in a state $s$ is $\rho_s$, the change of this
-+
+probability with time is given by
-+
+\begin{equation}
-+
+\frac{d \rho_s}{dt} = \sum_{s'} [ -w_{ss'}\rho_s + w_{s's}\rho_{s'} ],
-+
+\label{Ineq:timeChangeofProb}
-+
+\end{equation}
-+
+where $w_{ss'}$ is the tansition probability of going from state $s$
-+
+to state $s'$. Since $\rho_s$ is independent of time at equilibrium,
-+
+\begin{equation}
-+
+\frac{d \rho_{s}^{equilibrium}}{dt} = 0,
-+
+\label{Ineq:equiProb}
-+
+\end{equation}
-+
+which means $\sum_{s'} [ -w_{ss'}\rho_s + w_{s's}\rho_{s'} ]$ is $0$
-+
+for all $s'$. So
-+
+\begin{equation}
-+
+\frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = \frac{w_{s's}}{w_{ss'}}.
-+
+\label{Ineq:relationshipofRhoandW}
-+
+\end{equation}
-+
+If
-+
+\begin{equation}
-+
+\frac{w_{s's}}{w_{ss'}} = e^{-(U_s - U_{s'}) / k_B T},
-+
+\label{Ineq:conditionforBoltzmannStatistics}
-+
+\end{equation}
-+
+then
-+
+\begin{equation}
-+
+\frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = e^{-(U_s - U_{s'}) / k_B T}.
-+
+\label{Ineq:satisfyofBoltzmannStatistics}
-+
+\end{equation}
-+
+Eq.~\ref{Ineq:satisfyofBoltzmannStatistics} implies that
-+
+$\rho^{equilibrium}$ satisfies Boltzmann statistics. An algorithm,
-+
+shows how Monte Carlo simulation generates a transition probability
-+
+governed by \ref{Ineq:conditionforBoltzmannStatistics}, is schemed as
-+
+\begin{enumerate}
-+
+\item\label{Initm:oldEnergy} Choose an particle randomly, calculate the energy.
-+
+\item\label{Initm:newEnergy} Make a random displacement for particle,
-+
+calculate the new energy.
-+
+  \begin{itemize}
-+
+     \item Keep the new configuration and return to step~\ref{Initm:oldEnergy} if energy
-+
+goes down.
-+
+     \item Pick a random number between $[0,1]$ if energy goes up.
-+
+        \begin{itemize}
-+
+           \item Keep the new configuration and return to
-+
+step~\ref{Initm:oldEnergy} if the random number smaller than
-+
+$e^{-(U_{new} - U_{old})} / k_B T$.
-+
+           \item Keep the old configuration and return to
-+
+step~\ref{Initm:oldEnergy} if the random number larger than
-+
+$e^{-(U_{new} - U_{old})} / k_B T$.
-+
+        \end{itemize}
-+
+  \end{itemize}
-+
+\item\label{Initm:accumulateAvg} Accumulate the average after it converges.
-+
+\end{enumerate}
-+
+It is important to notice that the old configuration has to be sampled
-+
+again if it is kept.
-+
-+
+\subsection{Molecular Dynamics\label{In:ssec:md}}
-+
+Although some of properites of the system can be calculated from the
-+
+ensemble average in Monte Carlo simulations, due to the nature of
-+
+lacking in the time dependence, it is impossible to gain information
-+
+of those dynamic properties from Monte Carlo simulations. Molecular
-+
+Dynamics is a measurement of the evolution of the positions and
-+
+momenta of the particles in the system. The evolution of the system
-+
+obeys laws of classical mechanics, in most situations, there is no
-+
+need for the count of the quantum effects. For a real experiment, the
-+
+instantaneous positions and momenta of the particles in the system are
-+
+neither important nor measurable, the observable quantities are
-+
+usually a average value for a finite time interval. These quantities
-+
+are expressed as a function of positions and momenta in Melecular
-+
+Dynamics simulations. Like the thermal temperature of the system is
-+
+defined by
-+
+\begin{equation}
-+
+\frac{1}{2} k_B T = \langle \frac{1}{2} m v_\alpha \rangle,
-+
+\label{Ineq:temperature}
-+
+\end{equation}
-+
+here $m$ is the mass of the particle and $v_\alpha$ is the $\alpha$
-+
+component of the velocity of the particle. The right side of
-+
+Eq.~\ref{Ineq:temperature} is the average kinetic energy of the
-+
+system. A simple Molecular Dynamics simulation scheme
-+
+is:~\cite{Frenkel1996}
-+
+\begin{enumerate}
-+
+\item\label{Initm:initialize} Assign the initial positions and momenta
-+
+for the particles in the system.
-+
+\item\label{Initm:calcForce} Calculate the forces.
-+
+\item\label{Initm:equationofMotion} Integrate the equation of motion.
-+
+  \begin{itemize}
-+
+     \item Return to step~\ref{Initm:calcForce} if the equillibrium is
-+
+not achieved.
-+
+     \item Go to step~\ref{Initm:calcAvg} if the equillibrium is
-+
+achieved.
-+
+  \end{itemize}
-+
+\item\label{Initm:calcAvg} Compute the quantities we are interested in.
-+
+\end{enumerate}
-+
+The initial positions of the particles are chosen as that there is no
-+
+overlap for the particles. The initial velocities at first are
-+
+distributed randomly to the particles, and then shifted to make the
-+
+momentum of the system $0$, at last scaled to the desired temperature
-+
+of the simulation according Eq.~\ref{Ineq:temperature}.
-+
-+
+The core of Molecular Dynamics simulations is step~\ref{Initm:calcForce}
-+
+and~\ref{Initm:equationofMotion}. The calculation of the forces are
-+
+often involved numerous effort, this is the most time consuming step
-+
+in the Molecular Dynamics scheme. The evaluation of the forces is
-+
+followed by
-+
+\begin{equation}
-+
+f(q) = - \frac{\partial U(q)}{\partial q},
-+
+\label{Ineq:force}
-+
+\end{equation}
-+
+$U(q)$ is the potential of the system. Once the forces computed, are
-+
+the positions and velocities updated by integrating Newton's equation
-+
+of motion,
-+
+\begin{equation}
-+
+f(q) = \frac{dp}{dt} = \frac{m dv}{dt}.
-+
+\label{Ineq:newton}
-+
+\end{equation}
-+
+Here is an example of integrating algorithms, Verlet algorithm, which
-+
+is one of the best algorithms to integrate Newton's equation of
-+
+motion. The Taylor expension of the position at time $t$ is
-+
+\begin{equation}
-+
+q(t+\Delta t)= q(t) + v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2 +
-+
+        \frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} +
-+
+        \mathcal{O}(\Delta t^4)
-+
+\label{Ineq:verletFuture}
-+
+\end{equation}
-+
+for a later time $t+\Delta t$, and
-+
+\begin{equation}
-+
+q(t-\Delta t)= q(t) - v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2 -
-+
+        \frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} +
-+
+        \mathcal{O}(\Delta t^4) ,
-+
+\label{Ineq:verletPrevious}
-+
+\end{equation}
-+
+for a previous time $t-\Delta t$. The summation of the
-+
+Eq.~\ref{Ineq:verletFuture} and~\ref{Ineq:verletPrevious} gives
-+
+\begin{equation}
-+
+q(t+\Delta t)+q(t-\Delta t) =
-+
+q(t) + \frac{f(t)}{m}\Delta t^2 + \mathcal{O}(\Delta t^4),
-+
+\label{Ineq:verletSum}
-+
+\end{equation}
-+
+so, the new position can be expressed as
-+
+\begin{equation}
-+
+q(t+\Delta t) \approx
-+
+q(t) - q(t-\Delta t) + \frac{f(t)}{m}\Delta t^2.
-+
+\label{Ineq:newPosition}
-+
+\end{equation}
-+
+The higher order of the $\Delta t$ is omitted.
-+
-+
+Numerous technics and tricks are applied to Molecular Dynamics
-+
+simulation to gain more efficiency or more accuracy. The simulation
-+
+engine used in this dissertation for the Molecular Dynamics
-+
+simulations is {\sc oopse}, more details of the algorithms and
-+
+technics can be found in~\cite{Meineke2005}.
-+
-+
+\subsection{Langevin Dynamics\label{In:ssec:ld}}
-+
+In many cases, the properites of the solvent in a system, like the
-+
+lipid-water system studied in this dissertation, are less important to
-+
+the researchers. However, the major computational expense is spent on
-+
+the solvent in the Molecular Dynamics simulations because of the large
-+
+number of the solvent molecules compared to that of solute molecules,
-+
+the ratio of the number of lipid molecules to the number of water
-+
+molecules is $1:25$ in our lipid-water system. The efficiency of the
-+
+Molecular Dynamics simulations is greatly reduced.
-+
-+
+As an extension of the Molecular Dynamics simulations, the Langevin
-+
+Dynamics seeks a way to avoid integrating equation of motion for
-+
+solvent particles without losing the Brownian properites of solute
-+
+particles. A common approximation is that the coupling of the solute
-+
+and solvent is expressed as a set of harmonic oscillators. So the
-+
+Hamiltonian of such a system is written as
-+
+\begin{equation}
-+
+H = \frac{p^2}{2m} + U(q) + H_B + \Delta U(q),
-+
+\label{Ineq:hamiltonianofCoupling}
-+
+\end{equation}
-+
+where $H_B$ is the Hamiltonian of the bath which equals to
-+
+\begin{equation}
-+
+H_B = \sum_{\alpha = 1}^{N} \left\{ \frac{p_\alpha^2}{2m_\alpha} +
-+
+\frac{1}{2} m_\alpha \omega_\alpha^2 q_\alpha^2\right\},
-+
+\label{Ineq:hamiltonianofBath}
-+
+\end{equation}
-+
+$\alpha$ is all the degree of freedoms of the bath, $\omega$ is the
-+
+bath frequency, and $\Delta U(q)$ is the bilinear coupling given by
-+
+\begin{equation}
-+
+\Delta U = -\sum_{\alpha = 1}^{N} g_\alpha q_\alpha q,
-+
+\label{Ineq:systemBathCoupling}
-+
+\end{equation}
-+
+where $g$ is the coupling constant. By solving the Hamilton's equation
-+
+of motion, the {\it Generalized Langevin Equation} for this system is
-+
+derived to
-+
+\begin{equation}
-+
+m \ddot q = -\frac{\partial W(q)}{\partial q} - \int_0^t \xi(t) \dot q(t-t')dt' + R(t),
-+
+\label{Ineq:gle}
-+
+\end{equation}
-+
+with mean force,
-+
+\begin{equation}
-+
+W(q) = U(q) - \sum_{\alpha = 1}^N \frac{g_\alpha^2}{2m_\alpha
-+
+\omega_\alpha^2}q^2,
-+
+\label{Ineq:meanForce}
-+
+\end{equation}
-+
+being only a dependence of coordinates of the solute particles, {\it
-+
+friction kernel},
-+
+\begin{equation}
-+
+\xi(t) = \sum_{\alpha = 1}^N \frac{-g_\alpha}{m_\alpha
-+
+\omega_\alpha} \cos(\omega_\alpha t),
-+
+\label{Ineq:xiforGLE}
-+
+\end{equation}
-+
+and the random force,
-+
+\begin{equation}
-+
+R(t) = \sum_{\alpha = 1}^N \left( g_\alpha q_\alpha(0)-\frac{g_\alpha}{m_\alpha
-+
+\omega_\alpha^2}q(0)\right) \cos(\omega_\alpha t) + \frac{\dot
-+
+q_\alpha(0)}{\omega_\alpha} \sin(\omega_\alpha t),
-+
+\label{Ineq:randomForceforGLE}
-+
+\end{equation}
-+
+as only a dependence of the initial conditions. The relationship of
-+
+friction kernel $\xi(t)$ and random force $R(t)$ is given by
-+
+\begin{equation}
-+
+\xi(t) = \frac{1}{k_B T} \langle R(t)R(0) \rangle
-+
+\label{Ineq:relationshipofXiandR}
-+
+\end{equation}
-+
+from their definitions. In Langevin limit, the friction is treated
-+
+static, which means
-+
+\begin{equation}
-+
+\xi(t) = 2 \xi_0 \delta(t).
-+
+\label{Ineq:xiofStaticFriction}
-+
+\end{equation}
-+
+After substitude $\xi(t)$ into Eq.~\ref{Ineq:gle} with
-+
+Eq.~\ref{Ineq:xiofStaticFriction}, {\it Langevin Equation} is extracted
-+
+to
-+
+\begin{equation}
-+
+m \ddot q = -\frac{\partial U(q)}{\partial q} - \xi \dot q(t) + R(t).
-+
+\label{Ineq:langevinEquation}
-+
+\end{equation}
-+
+The applying of Langevin Equation to dynamic simulations is discussed
-+
+in Ch.~\ref{chap:ld}.

Diff Legend

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

Comparing trunk/xDissertation/Introduction.tex (file contents): Revision 3336 by xsun, Wed Jan 30 16:01:02 2008 UTC vs. Revision 3372 by xsun, Wed Mar 19 21:04:47 2008 UTC

Diff Legend

Comparing trunk/xDissertation/Introduction.tex (file contents):
Revision 3336 by xsun, Wed Jan 30 16:01:02 2008 UTC vs.
Revision 3372 by xsun, Wed Mar 19 21:04:47 2008 UTC