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 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}
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 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) = 
        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 (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}.
Revision:	3383
Committed:	Wed Apr 16 21:56:34 2008 UTC (17 years ago) by xsun
Content type:	application/x-tex
File size:	35348 byte(s)
Log Message:	formatting check