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. Lipids, when dispersed in water, self assemble into a
mumber of topologically distinct bilayer structures. The phase
behavior of lipid bilayers has been explored
experimentally~\cite{Cevc87}, however, a complete understanding of the
mechanism and driving forces behind the various phases has not been
achieved.

\subsection{Ripple Phase\label{In:ssec:ripple}}
The $P_{\beta'}$ {\it ripple phase} of lipid bilayers, named from the
periodic buckling of the membrane, is an intermediate phase which is
developed either from heating the gel phase $L_{\beta'}$ or cooling
the fluid phase $L_\alpha$. A Sketch is shown in
figure~\ref{Infig:phaseDiagram}.~\cite{Cevc87}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/inPhaseDiagram.pdf}
\caption{A phase diagram of lipid bilayer. With increasing the
temperature, the bilayer can go through a gel ($L_{\beta'}$), ripple
($P_{\beta'}$) to fluid ($L_\alpha$) phase transition.}
\label{Infig:phaseDiagram}
\end{figure}
Most structural information of the ripple phase has been obtained by
the X-ray diffraction~\cite{Sun96,Katsaras00} and freeze-fracture
electron microscopy (FFEM).~\cite{Copeland80,Meyer96} The X-ray
diffraction work by Katsaras {\it et al.} showed that a rich phase
diagram exhibiting both {\it asymmetric} and {\it symmetric} ripples
is possible for lecithin bilayers.\cite{Katsaras00} Recently,
Kaasgaard {\it et al.} used atomic force microscopy (AFM) to observe
ripple phase morphology in bilayers supported on
mica.~\cite{Kaasgaard03}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/inRipple.pdf}
\caption{The experimental observed ripple phase. The top image is
obtained by X-ray diffraction~\cite{Sun96}, and the bottom one is
observed by AFM.~\cite{Kaasgaard03}}
\label{Infig:ripple}
\end{figure}
Figure~\ref{Infig:ripple} shows the ripple phase oberved by X-ray
diffraction and AFM. The experimental results provide strong support
for a 2-dimensional triangular packing lattice of the lipid molecules
within the ripple phase.  This is a notable change from the observed
lipid packing within the gel phase,~\cite{Cevc87} although Tenchov
{\it et al.} have recently observed near-hexagonal packing in some
phosphatidylcholine (PC) gel phases.~\cite{Tenchov2001} However, the
physical mechanism for the formation of the ripple phase has never
been explained and the microscopic structure of the ripple phase has
never been elucidated by experiments. Computational simulation is a
perfect tool to study the microscopic properties for a
system. However, the large length scale the ripple structure and the
long time scale of the formation of the ripples are crucial obstacles
to performing the actual work. The principal ideas explored in this
dissertation are attempts to break the computational task up by
\begin{itemize}
\item Simplifying the lipid model.
\item Improving algorithm for integrating the equations of motion.
\end{itemize}
In chapters~\ref{chap:mc} and~\ref{chap:md}, we use a simple point
dipole spin model and a coarse-grained molecualr scale model to
perform the Monte Carlo and Molecular Dynamics simulations
respectively, and in chapter~\ref{chap:ld}, we develop a Langevin
Dynamics algorithm which excludes the explicit solvent to improve the
efficiency of the simulations.

\subsection{Lattice Model\label{In:ssec:model}}
The gel-like characteristic (relatively immobile molecules) exhibited
by the ripple phase makes it feasible to apply a lattice model to
study the system. The popular $2$ dimensional lattice models, {\it
i.e.}, the Ising, $X-Y$, and Heisenberg models, show {\it frustration}
on triangular lattices. The Hamiltonians of these systems are given by
\begin{equation}
H =
  \begin{cases}
    -J \sum_n \sum_{n'} s_n s_n' & \text{Ising}, \\
    -J \sum_n \sum_{n'} \vec s_n \cdot \vec s_{n'} & \text{$X-Y$ and
Heisenberg},
  \end{cases}
\end{equation}
where $J$ has non zero value only when spins $s_n$ ($\vec s_n$) and
$s_{n'}$ ($\vec s_{n'}$) are nearest neighbors. When $J > 0$, spins
prefer an aligned structure, and if $J < 0$, spins prefer an
anti-aligned structure.

\begin{figure}
\centering
\includegraphics[width=3in]{./figures/inFrustration.pdf}
\caption{Frustration on triangular lattice, the spins and dipoles are
represented by arrows. The multiple local minima of energy states
induce the frustration for spins and dipoles picking the directions.}
\label{Infig:frustration}
\end{figure}
The spins in figure~\ref{Infig:frustration} shows an illustration of
the frustration for $J < 0$ on a triangular lattice. There are
multiple local minima energy states which are independent of the
direction of the spin on top of the triangle, therefore infinite
possibilities for the packing of spins which induces what is known as
``complete regular frustration'' which leads to disordered low
temperature phases. The similarity goes to the dipoles on a hexagonal
lattice, which are shown by the dipoles in
figure~\ref{Infig:frustration}. In this circumstance, the dipoles want
to be aligned, however, due to the long wave fluctuation, at low
temperature, the aligned state becomes unstable, vortex is formed and
results in multiple local minima of energy states. The dipole on the
center of the hexagonal lattice is frustrated.

The lack of translational degrees of freedom in lattice models
prevents their utilization in models for surface buckling. In
chapter~\ref{chap:mc}, a modified lattice model is introduced to
tackle this specific situation.

\section{Overview of Classical Statistical Mechanics\label{In:sec:SM}}
Statistical mechanics provides a way to calculate the macroscopic
properties of a system from the molecular interactions used in
computational simulations. This section serves as a brief introduction
to key concepts of classical statistical mechanics that we used in
this dissertation. Tolman gives an excellent introduction to the
principles of statistical mechanics.~\cite{Tolman1979} A large part of
section~\ref{In:sec:SM} will follow Tolman's notation.

\subsection{Ensembles\label{In:ssec:ensemble}}
In classical mechanics, the state of the system is completely
described by the positions and momenta of all particles. If we have an
$N$ particle system, there are $6N$ coordinates ($3N$ positions $(q_1,
q_2, \ldots, q_{3N})$ and $3N$ momenta $(p_1, p_2, \ldots, p_{3N})$)
to define the instantaneous state of the system. Each single set of
the $6N$ coordinates can be considered as a unique point in a $6N$
dimensional space where each perpendicular axis is one of
$q_{i\alpha}$ or $p_{i\alpha}$ ($i$ is the particle and $\alpha$ is
the spatial axis). This $6N$ dimensional space is known as phase
space. The instantaneous state of the system is a single point in
phase space. A trajectory is a connected path of points in phase
space. An ensemble is a collection of systems described by the same
macroscopic observables but which have microscopic state
distributions. In phase space an ensemble is a collection of a set of
representive points. A density distribution $\rho(q^N, p^N)$ of the
representive points in phase space describes the condition of an
ensemble of identical systems. Since this density may change with
time, it is also a function of time. $\rho(q^N, p^N, t)$ describes the
ensemble at a time $t$, and $\rho(q^N, p^N, t')$ describes the
ensemble at a later time $t'$. For convenience, we will use $\rho$
instead of $\rho(q^N, p^N, t)$ in the following disccusion. If we
normalize $\rho$ to unity,
\begin{equation}
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), q^N(t+\tau)] \rho(q, p)
\label{Ineq:autocorrelationFunctionCM}
\end{equation}
and
\begin{equation}
C_{AB}(\tau) = \int d \vec q~^N \int d \vec p~^N A[(q^N(t), p^N(t)]
B[q^N(t+\tau), q^N(t+\tau)] \rho(q, p)
\label{Ineq:crosscorrelationFunctionCM}
\end{equation}
as autocorrelation function and cross correlation function
respectively. $\rho(q, p)$ is the density distribution at equillibrium
in phase space. In many cases, the correlation function decay is a
single exponential
\begin{equation}
C(t) \sim e^{-t / \tau_r},
\label{Ineq:relaxation}
\end{equation}
where $\tau_r$ is known as relaxation time which discribes the rate of
the decay.

\section{Methodolody\label{In:sec:method}}
The simulations performed in this dissertation are branched into two
main catalog, Monte Carlo and Molecular Dynamics. There are two main
difference between Monte Carlo and Molecular Dynamics simulations. One
is that the Monte Carlo simulation is time independent, and Molecular
Dynamics simulation is time involved. Another dissimilar is that the
Monte Carlo is a stochastic process, the configuration of the system
is not determinated by its past, however, using Moleuclar Dynamics,
the system is propagated by Newton's equation of motion, the
trajectory of the system evolved in the phase space is determined. A
brief introduction of the two algorithms are given in
section~\ref{In:ssec:mc} and~\ref{In:ssec:md}. An extension of the
Molecular Dynamics, Langevin Dynamics, is introduced by
section~\ref{In:ssec:ld}.

\subsection{Monte Carlo\label{In:ssec:mc}}
Monte Carlo algorithm was first introduced by Metropolis {\it et
al.}.~\cite{Metropolis53} Basic Monte Carlo algorithm is usually
applied to the canonical ensemble, a Boltzmann-weighted ensemble, in
which the $N$, the total number of particles, $V$, total volume, $T$,
temperature are constants. The average energy is given by substituding
Eq.~\ref{Ineq:canonicalEnsemble} into Eq.~\ref{Ineq:statAverage2},
\begin{equation}
\langle E \rangle = \frac{1}{Z_N} \int d \vec q~^N \int d \vec p~^N E e^{-H(q^N, p^N) / k_B T}.
\label{Ineq:energyofCanonicalEnsemble}
\end{equation}
So are the other properties of the system. The Hamiltonian is the
summation of Kinetic energy $K(p^N)$ as a function of momenta and
Potential energy $U(q^N)$ as a function of positions,
\begin{equation}
H(q^N, p^N) = K(p^N) + U(q^N).
\label{Ineq:hamiltonian}
\end{equation}
If the property $A$ is only a function of position ($ A = A(q^N)$),
the mean value of $A$ is reduced to
\begin{equation}
\langle A \rangle = \frac{\int d \vec q~^N \int d \vec p~^N A e^{-U(q^N) / k_B T}}{\int d \vec q~^N \int d \vec p~^N e^{-U(q^N) / k_B T}},
\label{Ineq:configurationIntegral}
\end{equation}
The kinetic energy $K(p^N)$ is factored out in
Eq.~\ref{Ineq:configurationIntegral}. $\langle A
\rangle$ is a configuration integral now, and the
Eq.~\ref{Ineq:configurationIntegral} is equivalent to
\begin{equation}
\langle A \rangle = \int d \vec q~^N A \rho(q^N).
\label{Ineq:configurationAve}
\end{equation}

In a Monte Carlo simulation of canonical ensemble, the probability of
the system being in a state $s$ is $\rho_s$, the change of this
probability with time is given by
\begin{equation}
\frac{d \rho_s}{dt} = \sum_{s'} [ -w_{ss'}\rho_s + w_{s's}\rho_{s'} ],
\label{Ineq:timeChangeofProb}
\end{equation}
where $w_{ss'}$ is the tansition probability of going from state $s$
to state $s'$. Since $\rho_s$ is independent of time at equilibrium,
\begin{equation}
\frac{d \rho_{s}^{equilibrium}}{dt} = 0,
\label{Ineq:equiProb}
\end{equation}
which means $\sum_{s'} [ -w_{ss'}\rho_s + w_{s's}\rho_{s'} ]$ is $0$
for all $s'$. So
\begin{equation}
\frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = \frac{w_{s's}}{w_{ss'}}.
\label{Ineq:relationshipofRhoandW}
\end{equation}
If
\begin{equation}
\frac{w_{s's}}{w_{ss'}} = e^{-(U_s - U_{s'}) / k_B T},
\label{Ineq:conditionforBoltzmannStatistics}
\end{equation}
then
\begin{equation}
\frac{\rho_s^{equilibrium}}{\rho_{s'}^{equilibrium}} = e^{-(U_s - U_{s'}) / k_B T}.
\label{Ineq:satisfyofBoltzmannStatistics}
\end{equation}
Eq.~\ref{Ineq:satisfyofBoltzmannStatistics} implies that
$\rho^{equilibrium}$ satisfies Boltzmann statistics. An algorithm,
shows how Monte Carlo simulation generates a transition probability
governed by \ref{Ineq:conditionforBoltzmannStatistics}, is schemed as
\begin{enumerate}
\item\label{Initm:oldEnergy} Choose an particle randomly, calculate the energy.
\item\label{Initm:newEnergy} Make a random displacement for particle,
calculate the new energy.
  \begin{itemize}
     \item Keep the new configuration and return to step~\ref{Initm:oldEnergy} if energy
goes down.
     \item Pick a random number between $[0,1]$ if energy goes up.
        \begin{itemize}
           \item Keep the new configuration and return to
step~\ref{Initm:oldEnergy} if the random number smaller than
$e^{-(U_{new} - U_{old})} / k_B T$.
           \item Keep the old configuration and return to
step~\ref{Initm:oldEnergy} if the random number larger than
$e^{-(U_{new} - U_{old})} / k_B T$.
        \end{itemize}
  \end{itemize}
\item\label{Initm:accumulateAvg} Accumulate the average after it converges.
\end{enumerate}
It is important to notice that the old configuration has to be sampled
again if it is kept.

\subsection{Molecular Dynamics\label{In:ssec:md}}
Although some of properites of the system can be calculated from the
ensemble average in Monte Carlo simulations, due to the nature of
lacking in the time dependence, it is impossible to gain information
of those dynamic properties from Monte Carlo simulations. Molecular
Dynamics is a measurement of the evolution of the positions and
momenta of the particles in the system. The evolution of the system
obeys laws of classical mechanics, in most situations, there is no
need for the count of the quantum effects. For a real experiment, the
instantaneous positions and momenta of the particles in the system are
neither important nor measurable, the observable quantities are
usually a average value for a finite time interval. These quantities
are expressed as a function of positions and momenta in Melecular
Dynamics simulations. Like the thermal temperature of the system is
defined by
\begin{equation}
\frac{1}{2} k_B T = \langle \frac{1}{2} m v_\alpha \rangle,
\label{Ineq:temperature}
\end{equation}
here $m$ is the mass of the particle and $v_\alpha$ is the $\alpha$
component of the velocity of the particle. The right side of
Eq.~\ref{Ineq:temperature} is the average kinetic energy of the
system. A simple Molecular Dynamics simulation scheme
is:~\cite{Frenkel1996}
\begin{enumerate}
\item\label{Initm:initialize} Assign the initial positions and momenta
for the particles in the system.
\item\label{Initm:calcForce} Calculate the forces.
\item\label{Initm:equationofMotion} Integrate the equation of motion.
  \begin{itemize}
     \item Return to step~\ref{Initm:calcForce} if the equillibrium is
not achieved.
     \item Go to step~\ref{Initm:calcAvg} if the equillibrium is
achieved.
  \end{itemize}
\item\label{Initm:calcAvg} Compute the quantities we are interested in.
\end{enumerate}
The initial positions of the particles are chosen as that there is no
overlap for the particles. The initial velocities at first are
distributed randomly to the particles, and then shifted to make the
momentum of the system $0$, at last scaled to the desired temperature
of the simulation according Eq.~\ref{Ineq:temperature}.

The core of Molecular Dynamics simulations is step~\ref{Initm:calcForce}
and~\ref{Initm:equationofMotion}. The calculation of the forces are
often involved numerous effort, this is the most time consuming step
in the Molecular Dynamics scheme. The evaluation of the forces is
followed by
\begin{equation}
f(q) = - \frac{\partial U(q)}{\partial q},
\label{Ineq:force}
\end{equation}
$U(q)$ is the potential of the system. Once the forces computed, are
the positions and velocities updated by integrating Newton's equation
of motion,
\begin{equation}
f(q) = \frac{dp}{dt} = \frac{m dv}{dt}.
\label{Ineq:newton}
\end{equation}
Here is an example of integrating algorithms, Verlet algorithm, which
is one of the best algorithms to integrate Newton's equation of
motion. The Taylor expension of the position at time $t$ is
\begin{equation}
q(t+\Delta t)= q(t) + v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2 +
        \frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} +
        \mathcal{O}(\Delta t^4)
\label{Ineq:verletFuture}
\end{equation}
for a later time $t+\Delta t$, and 
\begin{equation}
q(t-\Delta t)= q(t) - v(t) \Delta t + \frac{f(t)}{2m}\Delta t^2 -
        \frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} +
        \mathcal{O}(\Delta t^4) ,
\label{Ineq:verletPrevious}
\end{equation}
for a previous time $t-\Delta t$. The summation of the
Eq.~\ref{Ineq:verletFuture} and~\ref{Ineq:verletPrevious} gives
\begin{equation}
q(t+\Delta t)+q(t-\Delta t) = 
        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 of the $\Delta t$ is omitted.

Numerous technics and tricks are applied to Molecular Dynamics
simulation to gain more efficiency or more accuracy. The simulation
engine used in this dissertation for the Molecular Dynamics
simulations is {\sc oopse}, more details of the algorithms and
technics can be found in~\cite{Meineke2005}.

\subsection{Langevin Dynamics\label{In:ssec:ld}}
In many cases, the properites of the solvent in a system, like the
lipid-water system studied in this dissertation, are less important to
the researchers. However, the major computational expense is spent on
the solvent in the Molecular Dynamics simulations because of the large
number of the solvent molecules compared to that of solute molecules,
the ratio of the number of lipid molecules to the number of water
molecules is $1:25$ in our lipid-water system. The efficiency of the
Molecular Dynamics simulations is greatly reduced.

As an extension of the Molecular Dynamics simulations, the Langevin
Dynamics seeks a way to avoid integrating equation of motion for
solvent particles without losing the Brownian properites of solute
particles. A common approximation is that the coupling of the solute
and solvent is expressed as a set of harmonic oscillators. So the
Hamiltonian of such a system is written as
\begin{equation}
H = \frac{p^2}{2m} + U(q) + H_B + \Delta U(q),
\label{Ineq:hamiltonianofCoupling}
\end{equation}
where $H_B$ is the Hamiltonian of the bath which equals to
\begin{equation}
H_B = \sum_{\alpha = 1}^{N} \left\{ \frac{p_\alpha^2}{2m_\alpha} +
\frac{1}{2} m_\alpha \omega_\alpha^2 q_\alpha^2\right\},
\label{Ineq:hamiltonianofBath}
\end{equation}
$\alpha$ is all the degree of freedoms of the bath, $\omega$ is the
bath frequency, and $\Delta U(q)$ is the bilinear coupling given by
\begin{equation}
\Delta U = -\sum_{\alpha = 1}^{N} g_\alpha q_\alpha q,
\label{Ineq:systemBathCoupling}
\end{equation}
where $g$ is the coupling constant. By solving the Hamilton's equation
of motion, the {\it Generalized Langevin Equation} for this system is
derived to
\begin{equation}
m \ddot q = -\frac{\partial W(q)}{\partial q} - \int_0^t \xi(t) \dot q(t-t')dt' + R(t),
\label{Ineq:gle}
\end{equation}
with mean force, 
\begin{equation}
W(q) = U(q) - \sum_{\alpha = 1}^N \frac{g_\alpha^2}{2m_\alpha
\omega_\alpha^2}q^2,
\label{Ineq:meanForce}
\end{equation}
being only a dependence of coordinates of the solute particles, {\it
friction kernel},
\begin{equation}
\xi(t) = \sum_{\alpha = 1}^N \frac{-g_\alpha}{m_\alpha
\omega_\alpha} \cos(\omega_\alpha t),
\label{Ineq:xiforGLE}
\end{equation}
and the random force,
\begin{equation}
R(t) = \sum_{\alpha = 1}^N \left( g_\alpha q_\alpha(0)-\frac{g_\alpha}{m_\alpha
\omega_\alpha^2}q(0)\right) \cos(\omega_\alpha t) + \frac{\dot
q_\alpha(0)}{\omega_\alpha} \sin(\omega_\alpha t),
\label{Ineq:randomForceforGLE}
\end{equation}
as only a dependence of the initial conditions. The relationship of
friction kernel $\xi(t)$ and random force $R(t)$ is given by
\begin{equation}
\xi(t) = \frac{1}{k_B T} \langle R(t)R(0) \rangle
\label{Ineq:relationshipofXiandR}
\end{equation}
from their definitions. In Langevin limit, the friction is treated
static, which means
\begin{equation}
\xi(t) = 2 \xi_0 \delta(t).
\label{Ineq:xiofStaticFriction}
\end{equation}
After substitude $\xi(t)$ into Eq.~\ref{Ineq:gle} with
Eq.~\ref{Ineq:xiofStaticFriction}, {\it Langevin Equation} is extracted
to
\begin{equation}
m \ddot q = -\frac{\partial U(q)}{\partial q} - \xi \dot q(t) + R(t).
\label{Ineq:langevinEquation}
\end{equation}
The applying of Langevin Equation to dynamic simulations is discussed
in Ch.~\ref{chap:ld}.
Revision:	3372
Committed:	Wed Mar 19 21:04:47 2008 UTC (17 years, 1 month ago) by xsun
Content type:	application/x-tex
File size:	32509 byte(s)
Log Message:	writing.