[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 3383 by xsun, Wed Apr 16 21:56:34 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. The chemical structure of phospholipids includes a head
-+
+group with a large dipole moment which is due to the large charge
-+
+separation between phosphate and amino alcohol, and a nonpolar tail
-+
+that contains fatty acid chains. Depending on the specific alcohol
-+
+which the phosphate and fatty acid chains are esterified to, the
-+
+phospholipids are divided into glycerophospholipids and
-+
+sphingophospholipids.~\cite{Cevc80} The chemical structures are shown
-+
+in figure~\ref{Infig:lipid}.
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=\linewidth]{./figures/inLipid.pdf}
-+
+\caption[The chemical structure of lipids]{The chemical structure of
-+
+glycerophospholipids (left) and sphingophospholipids
-+
+(right).\cite{Cevc80}}
-+
+\label{Infig:lipid}
-+
+\end{figure}
-+
+Glycerophospholipids are the dominant phospholipids in biological
-+
+membranes. The type of glycerophospholipid depends on the identity of
-+
+the X group, and the chains. For example, if X is choline
-+
+[(CH$_3$)$_3$N$^+$CH$_2$CH$_2$OH], the lipids are known as
-+
+phosphatidylcholine (PC), or if X is ethanolamine
-+
+[H$_3$N$^+$CH$_2$CH$_2$OH], the lipids are known as
-+
+phosphatidyethanolamine (PE). Table~\ref{Intab:pc} listed a number
-+
+types of phosphatidycholine with different fatty acids as the lipid
-+
+chains.
-+
+\begin{table*}
-+
+\begin{minipage}{\linewidth}
-+
+\begin{center}
-+
+\caption{A NUMBER TYPES OF PHOSPHATIDYCHOLINE}
-+
+\begin{tabular}{lll}
-+
+\hline
-+
+  & Fatty acid & Generic Name \\
-+
+\hline
-+
+\textcolor{red}{DMPC} & Myristic: CH$_3$(CH$_2$)$_{12}$COOH &
-+
+\textcolor{red}{D}i\textcolor{red}{M}yristoyl\textcolor{red}{P}hosphatidyl\textcolor{red}{C}holine \\
-+
+\textcolor{red}{DPPC} & Palmitic: CH$_3$(CH$_2$)$_{14}$COOH & \textcolor{red}{D}i\textcolor{red}{P}almtoyl\textcolor{red}{P}hosphatidyl\textcolor{red}{C}holine
-+
+\\
-+
+\textcolor{red}{DSPC} & Stearic: CH$_3$(CH$_2$)$_{16}$COOH & \textcolor{red}{D}i\textcolor{red}{S}tearoyl\textcolor{red}{P}hosphatidyl\textcolor{red}{C}holine \\
-+
+\end{tabular}
-+
+\label{Intab:pc}
-+
+\end{center}
-+
+\end{minipage}
-+
+\end{table*}
-+
+When dispersed in water, lipids self assemble into a number of
-+
+topologically distinct bilayer structures. The phase behavior of lipid
-+
+bilayers has been explored experimentally~\cite{Cevc80}, however, a
-+
+complete understanding of the mechanism and driving forces behind the
-+
+various phases has not been achieved.
-+
-+
+\subsection{The 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 of the phases is shown in
-+
+figure~\ref{Infig:phaseDiagram}.~\cite{Cevc80}
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=\linewidth]{./figures/inPhaseDiagram.pdf}
-+
+\caption[Phases of PC lipid bilayers]{Phases of PC lipid
-+
+bilayers. With increasing temperature, phosphotidylcholine (PC)
-+
+bilayers can go through $L_{\beta'} \rightarrow P_{\beta'}$ (gel
-+
+$\rightarrow$ ripple) and $P_{\beta'} \rightarrow L_\alpha$ (ripple
-+
+$\rightarrow$ fluid) phase transitions.~\cite{Cevc80}}
-+
+\label{Infig:phaseDiagram}
-+
+\end{figure}
-+
+Most structural information about the ripple phase has been obtained
-+
+by 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[Experimental observations of the riple phase]{Experimental
-+
+observations of the riple phase. The top image is an electrostatic
-+
+density map obtained by Sun {\it et al.} using X-ray
-+
+diffraction~\cite{Sun96}.  The lower figures are the surface topology
-+
+of various ripple domains in bilayers supported in mica. The AFM
-+
+images were observed by Kaasgaard {\it et al.}.~\cite{Kaasgaard03}}
-+
+\label{Infig:ripple}
-+
+\end{figure}
-+
+Figure~\ref{Infig:ripple} shows the ripple phase oberved by both 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{Cevc80} although Tenchov
-+
+{\it et al.} have recently observed near-triangular 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
-+
+very good tool to study the microscopic properties for a
-+
+system. However, the large length scale of the ripple structures and
-+
+the long time required for 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 the 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 molecular 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 Models\label{In:ssec:model}}
-+
+The gel-like characteristic (laterally 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]{Frustration on triangular
-+
+lattice, the spins and dipoles are represented by arrows. The multiple
-+
+local minima of energy states induce frustration for spins and dipoles
-+
+resulting in disordered low-temperature phases.}
-+
+\label{Infig:frustration}
-+
+\end{figure}
-+
+The spins in figure~\ref{Infig:frustration} illustrate 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 orienting
-+
+large numbers spins. This induces what is known as ``complete regular
-+
+frustration'' which leads to disordered low temperature phases. This
-+
+behavior extends to dipoles on a triangular lattices, which are shown
-+
+in the lower portion of figure~\ref{Infig:frustration}. In this case,
-+
+dipole-aligned structures are energetically favorable, however, at low
-+
+temperature, vortices are easily formed, and, this results in multiple
-+
+local minima of energy states for a central dipole. The dipole on the
-+
+center of the hexagonal lattice is therefore 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} follows 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 type of thermodynamic 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}, but this ensemble has
-+
+little relevance for physical chemistry. It is an ensemble with
-+
+essentially infinite temperature.
-+
-+
+\subsubsection{The Microcanonical Ensemble\label{In:sssec:microcanonical}}
-+
+The most useful ensemble for Molecular Dynamics 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 which has the same units as
-+
+$\Sigma(N, V, E)$. 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. 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. Tow representative points in phase space are equally likely to be
-+
+found if they have the same energy. In other words, all energetically
-+
+accessible states are represented , and the probabilities to find the
-+
+system in any of the accessible states is equal to that states
-+
+Boltzmann weight. 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 finite time interval, $T$.
-+
-+
+\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 is called a {\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 mechanical form,
-+
+\begin{equation}
-+
+C_{AA}(\tau) = \int d \vec q~^N \int d \vec p~^N A[(q^N, p^N]
-+
+A[q^N(\tau), p^N(\tau)] \rho(q^N, p^N)
-+
+\label{Ineq:autocorrelationFunctionCM}
-+
+\end{equation}
-+
+where $q^N(\tau)$, $p^N(\tau)$ is the phase space point that follows
-+
+classical evolution of the point $q^N$, $p^N$ after a tme $\tau$ has
-+
+elapsed, and
-+
+\begin{equation}
-+
+C_{AB}(\tau) = \int d \vec q~^N \int d \vec p~^N A[(q^N, p^N]
-+
+B[q^N(\tau), p^N(\tau)] \rho(q^N, p^N)
-+
+\label{Ineq:crosscorrelationFunctionCM}
-+
+\end{equation}
-+
+as the autocorrelation function and cross correlation functions
-+
+respectively. $\rho(q, p)$ is the density distribution at equillibrium
-+
+in phase space. In many cases, correlation functions decay as a
-+
+single exponential in time
-+
+\begin{equation}
-+
+C(t) \sim e^{-t / \tau_r},
-+
+\label{Ineq:relaxation}
-+
+\end{equation}
-+
+where $\tau_r$ is known as relaxation time which describes the rate of
-+
+the decay.
-+
-+
+\section{Methodology\label{In:sec:method}}
-+
+The simulations performed in this dissertation branch into two main
-+
+categories, Monte Carlo and Molecular Dynamics. There are two main
-+
+differences between Monte Carlo and Molecular Dynamics
-+
+simulations. One is that the Monte Carlo simulations are time
-+
+independent methods of sampling structural features of an ensemble,
-+
+while Molecular Dynamics simulations provide dynamic
-+
+information. Additionally, Monte Carlo methods are stochastic
-+
+processes; the future configurations of the system are not determined
-+
+by its past. However, in Molecular Dynamics, the system is propagated
-+
+by Hamilton's equations of motion, and the trajectory of the system
-+
+evolving in phase space is deterministic. Brief introductions of the
-+
+two algorithms are given in section~\ref{In:ssec:mc}
-+
+and~\ref{In:ssec:md}. Langevin Dynamics, an extension of the Molecular
-+
+Dynamics that includes implicit solvent effects, is introduced by
-+
+section~\ref{In:ssec:ld}.
-+
-+
+\subsection{Monte Carlo\label{In:ssec:mc}}
-+
+A Monte Carlo integration algorithm was first introduced by Metropolis
-+
+{\it et al.}~\cite{Metropolis53} The basic Metropolis Monte Carlo
-+
+algorithm is usually applied to the canonical ensemble, a
-+
+Boltzmann-weighted ensemble, in which $N$, the total number of
-+
+particles, $V$, the total volume, and $T$, the temperature are
-+
+constants. An average in this ensemble is given
-+
+\begin{equation}
-+
+\langle A \rangle = \frac{1}{Z_N} \int d \vec q~^N \int d \vec p~^N
-+
+A(q^N, p^N) e^{-H(q^N, p^N) / k_B T}.
-+
+\label{Ineq:energyofCanonicalEnsemble}
-+
+\end{equation}
-+
+The Hamiltonian is the sum of the kinetic energy, $K(p^N)$, a function
-+
+of momenta and the potential energy, $U(q^N)$, 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 a function only of position ($ A = A(q^N)$),
-+
+the mean value of $A$ can be reduced to
-+
+\begin{equation}
-+
+\langle A \rangle = \frac{\int d \vec q~^N A e^{-U(q^N) / k_B T}}{\int d \vec q~^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 now a configuration integral, and
-+
+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}
-+
+where $\rho(q^N)$ is a configurational probability
-+
+\begin{equation}
-+
+\rho(q^N) = \frac{e^{-U(q^N) / k_B T}}{\int d \vec q~^N e^{-U(q^N) / k_B T}}.
-+
+\label{Ineq:configurationProb}
-+
+\end{equation}
-+
-+
+In a Monte Carlo simulation of a system in the 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}
-+
+the sum of transition probabilities $\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 the ratio of state populations
-+
+\begin{equation}
-+
+\frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = e^{-(U_s - U_{s'}) / k_B T}.
-+
+\label{Ineq:satisfyofBoltzmannStatistics}
-+
+\end{equation}
-+
+then the ratio of transition probabilities,
-+
+\begin{equation}
-+
+\frac{w_{s's}}{w_{ss'}} = e^{-(U_s - U_{s'}) / k_B T},
-+
+\label{Ineq:conditionforBoltzmannStatistics}
-+
+\end{equation}
-+
+An algorithm that indicates how a Monte Carlo simulation generates a
-+
+transition probability governed by
-+
+\ref{Ineq:conditionforBoltzmannStatistics}, is given schematically as,
-+
+\begin{enumerate}
-+
+\item\label{Initm:oldEnergy} Choose a particle randomly, and calculate
-+
+the energy of the rest of the system due to the current location of
-+
+the particle.
-+
+\item\label{Initm:newEnergy} Make a random displacement of the particle,
-+
+calculate the new energy taking the movement of the particle into account.
-+
+  \begin{itemize}
-+
+     \item If the energy goes down, keep the new configuration.
-+
+     \item If the energy goes up, pick a random number between $[0,1]$.
-+
+        \begin{itemize}
-+
+           \item If the random number smaller than
-+
+$e^{-(U_{new} - U_{old})} / k_B T$, keep the new configuration.
-+
+           \item If the random number is larger than
-+
+$e^{-(U_{new} - U_{old})} / k_B T$, keep the old configuration.
-+
+        \end{itemize}
-+
+  \end{itemize}
-+
+\item\label{Initm:accumulateAvg} Accumulate the averages based on the
-+
+current configuration.
-+
+\item Go to step~\ref{Initm:oldEnergy}.
-+
+\end{enumerate}
-+
+It is important for sampling accuracy that the old configuration is
-+
+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 absence of the
-+
+time dependence, it is impossible to gain information on dynamic
-+
+properties from Monte Carlo simulations. Molecular Dynamics evolves
-+
+the positions and momenta of the particles in the system. The
-+
+evolution of the system obeys the laws of classical mechanics, and in
-+
+most situations, there is no need to account for quantum effects. In a
-+
+real experiment, the instantaneous positions and momenta of the
-+
+particles in the system are ofter neither important nor measurable,
-+
+the observable quantities are usually an average value for a finite
-+
+time interval. These quantities are expressed as a function of
-+
+positions and momenta in Molecular Dynamics simulations. For example,
-+
+temperature of the system is defined by
-+
+\begin{equation}
-+
+\frac{3}{2} N k_B T = \langle \sum_{i=1}^N \frac{1}{2} m_i v_i \rangle,
-+
+\label{Ineq:temperature}
-+
+\end{equation}
-+
+here $m_i$ is the mass of particle $i$ and $v_i$ is the velocity of
-+
+particle $i$. The right side of Eq.~\ref{Ineq:temperature} is the
-+
+average kinetic energy of the system.
-+
-+
+The initial positions of the particles are chosen so that there is no
-+
+overlap of the particles. The initial velocities at first are
-+
+distributed randomly to the particles using a Maxwell-Boltzmann
-+
+distribution, and then shifted to make the total linear momentum of
-+
+the system $0$.
-+
-+
+The core of Molecular Dynamics simulations is the calculation of
-+
+forces and the integration algorithm. Calculation of the forces often
-+
+involves enormous effort. This is the most time consuming step in the
-+
+Molecular Dynamics scheme. Evaluation of the forces is mathematically
-+
+simple,
-+
+\begin{equation}
-+
+f(q) = - \frac{\partial U(q)}{\partial q},
-+
+\label{Ineq:force}
-+
+\end{equation}
-+
+where $U(q)$ is the potential of the system. However, the numerical
-+
+details of this computation are often quite complex. Once the forces
-+
+computed, the positions and velocities are updated by integrating
-+
+Hamilton's equations of motion,
-+
+\begin{eqnarray}
-+
+\dot p_i & = & -\frac{\partial H}{\partial q_i} = -\frac{\partial
-+
+U(q_i)}{\partial q_i} = f(q_i) \\
-+
+\dot q_i & = & p_i
-+
+\label{Ineq:newton}
-+
+\end{eqnarray}
-+
+The classic example of an integrating algorithm is the Verlet
-+
+algorithm, which is one of the simplest algorithms for integrating the
-+
+equations of motion. The Taylor expansion of the position at time $t$
-+
+is
-+
+\begin{equation}
-+
+q(t+\Delta t)= q(t) + \frac{p(t)}{m} \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) - \frac{p(t)}{m} \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$. Adding 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 terms in $\Delta t$ are omitted.
-+
-+
+Numerous techniques and tricks have been applied to Molecular Dynamics
-+
+simulations to gain greater efficiency or accuracy. The engine used in
-+
+this dissertation for the Molecular Dynamics simulations is {\sc
-+
+oopse}. More details of the algorithms and techniques used in this
-+
+code can be found in Ref.~\cite{Meineke2005}.
-+
-+
+\subsection{Langevin Dynamics\label{In:ssec:ld}}
-+
+In many cases, the properites of the solvent (like the water in the
-+
+lipid-water system studied in this dissertation) are less interesting
-+
+to the researchers than the behavior of the solute. However, the major
-+
+computational expense is ofter the solvent-solvent interactions, this
-+
+is due to the large number of the solvent molecules when compared to
-+
+the number 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, with up to 85\% of CPU time spent
-+
+calculating only water-water interactions.
-+
-+
+As an extension of the Molecular Dynamics methodologies, Langevin
-+
+Dynamics seeks a way to avoid integrating the equations of motion for
-+
+solvent particles without losing the solvent effects on the solute
-+
+particles. One common approximation is to express the coupling of the
-+
+solute and solvent as a set of harmonic oscillators. 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  is a set of N
-+
+harmonic oscillators
-+
+\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$ runs over all the degree of freedoms of the bath,
-+
+$\omega_\alpha$ is the bath frequency of oscillator $\alpha$, 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_\alpha$ is the coupling constant for oscillator $\alpha$. By
-+
+solving the Hamilton's equations of motion, the {\it Generalized
-+
+Langevin Equation} for this system is derived as
-+
+\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}
-+
+The {\it friction kernel}, $\xi(t)$, depends only on the coordinates
-+
+of the solute particles,
-+
+\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 a ``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}
-+
+that depends only on the initial conditions. The relationship of
-+
+friction kernel $\xi(t)$ and random force $R(t)$ is given by the
-+
+second fluctuation dissipation theorem,
-+
+\begin{equation}
-+
+\xi(t) = \frac{1}{k_B T} \langle R(t)R(0) \rangle.
-+
+\label{Ineq:relationshipofXiandR}
-+
+\end{equation}
-+
+In the harmonic bath this relation is exact and provable from the
-+
+definitions of these quantities. In the limit of static friction,
-+
+\begin{equation}
-+
+\xi(t) = 2 \xi_0 \delta(t).
-+
+\label{Ineq:xiofStaticFriction}
-+
+\end{equation}
-+
+After substituting $\xi(t)$ into Eq.~\ref{Ineq:gle} with
-+
+Eq.~\ref{Ineq:xiofStaticFriction}, the {\it Langevin Equation} is
-+
+extracted,
-+
+\begin{equation}
-+
+m \ddot q = -\frac{\partial U(q)}{\partial q} - \xi \dot q(t) + R(t).
-+
+\label{Ineq:langevinEquation}
-+
+\end{equation}
-+
+Application of the 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 3383 by xsun, Wed Apr 16 21:56:34 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 3383 by xsun, Wed Apr 16 21:56:34 2008 UTC