trunk/xDissertation/Introduction.tex

\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 the polar
head group which is due to a large charge separation between phosphate
and amino alcohol, and the nonpolar tails that contains fatty acid
chains. Depending on the alcohol which 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 glycerophospholipids (left) and
sphingophospholipids (right).\cite{Cevc80}}
\label{Infig:lipid}
\end{figure}
The glycerophospholipid is the dominant phospholipid in membranes. The
types of glycerophospholipids depend on the X group, and the
chains. For example, if X is choline, the lipids are known as
phosphatidylcholine (PC), or if X is ethanolamine, 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 mumber 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{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{Cevc80}
\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.~\cite{Cevc80}}
\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 Sun {\it et al.} using X-ray diffraction~\cite{Sun96},
and the bottom one is observed by Kaasgaard {\it et al.} using
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{Cevc80} 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}
1 = \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), p^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), p^N(t+\tau)] \rho(q, p)
\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 a stochastic
processes, the future configurations of the system are not determined
by its past. However, in Molecular Dynamics, the system is propagated
by Newton's equation 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 the $N$, the total number of
particles, $V$, total volume, $T$, 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}

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}
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 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 shows how 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 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 and return to
step~\ref{Initm:oldEnergy}.
     \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 and return to
step~\ref{Initm:oldEnergy}.
           \item If the random number is larger than
$e^{-(U_{new} - U_{old})} / k_B T$, keep the old configuration and return to
step~\ref{Initm:oldEnergy}.
        \end{itemize}
  \end{itemize}
\item\label{Initm:accumulateAvg} Accumulate the averages based on the
current configuartion.
\item Go to step~\ref{Initm:oldEnergy}.
\end{enumerate}
It is important for sampling accurately 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 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
ditribution, 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) + 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$. Adding Eq.~\ref{Ineq:verletFuture}
and~\ref{Ineq:verletPrevious} gives
\begin{equation}
q(t+\Delta t)+q(t-\Delta t) = 
        2q(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 
        2q(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 in a system, 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 interation, 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}
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}.
Revision:	3374
Committed:	Thu Mar 20 22:21:43 2008 UTC (17 years, 7 months ago) by xsun
Content type:	application/x-tex
File size:	34376 byte(s)
Log Message:	writing.