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 bilayer
structures. The phase behavior of lipid bilayers have been explored
experimentally~\cite{Cevc87}, however, a complete understanding of the
mechanism for self-assembly 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$. Although the ripple phase has been
observed in several
experiments~\cite{Sun96,Katsaras00,Copeland80,Meyer96,Kaasgaard03},
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 model and a coarse-grained 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 of the ripple phase makes it feasible to
apply a lattice model to study the system. It is claimed that the
packing of the lipid molecules in ripple phase is
hexagonal~\cite{Cevc87}. The popular $2$ dimensional lattice models,
{\it i.e.}, the Ising, $X-Y$, and Heisenberg models, show {\it
frustration} on triangular lattices. The Hamiltonian of the 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 the nearest neighbours.
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/inFrustration.pdf}
\caption{Frustration on a triangular lattice, the spins are
represented by arrows. No matter which direction the spin on the top
of triangle points to, the Hamiltonain of the system goes up.}
\label{Infig:frustration}
\end{figure}
Figure~\ref{Infig:frustration} shows an illustration of the
frustration on a triangular lattice. When $J < 0$, the spins want to
be anti-aligned, The direction of the spin on top of the triangle will
make the energy go up no matter which direction it picks, therefore
infinite possibilities for the packing of spins induce what is known
as ``complete regular frustration'' which leads to disordered low
temperature phases.

The lack of translational degree of freedom in lattice models prevents
their utilization in models for surface buckling which would
correspond to ripple formation. 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$ position $(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 the phase space.

Liouville's theorem describes the change in density $\rho$ with
time. The number of representive points at a given volume in the 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 representive points in this
volume, let us consider a simple condition: the change in the number
of representive points in $q_1$ axis. The rate of the number of the
representive 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 representive 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 of
the number of the representive points is the difference of
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 of the number of
representive 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 equation 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
can derive 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 of $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 representive 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 representive 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 often 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. Here are several
examples: 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}.  Another useful
ensemble is called {\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, it is the
total phase space volume of one state. The entropy in 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}
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 {\it partition function}. $\Gamma$ indicates
that the integral is over all the 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,
the thermodynamics maximizes the entropy $S$,
\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 of the canonical
ensemble, {\it i.e.}, total probability and average energy conserved,
the partition function is calculated as
\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 the microscopic states.

There is an implicit assumption that our arguments are based on so
far. A representive point in the phase space is equally to be found in
any same extent of the regions. 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
equal a {\it priori} probabilities.

\subsection{Statistical Average\label{In:ssec:average}}
Given the density distribution $\rho$ in the 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, t) \rangle = \frac{\int d \vec q~^N \int d \vec p~^N
F(q^N, p^N, t) \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, t) \rangle = \int d \vec q~^N \int d \vec p~^N F(q^N,
p^N, t) \rho,
\label{Ineq:statAverage2}
\end{equation}
called {\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, 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, p^N, t) dt
\label{Ineq:timeAverage2}
\end{equation}
for an infinite time interval.

{\it ergodic hypothesis}, an important hypothesis from the statistical
mechanics point of view, 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, t) \rangle = \langle F(q^N, p^N) \rangle_t.
\label{Ineq:ergodicity}
\end{equation}
Eq.~\ref{Ineq:ergodicity} validates the Monte Carlo method which we will
discuss in section~\ref{In:ssec:mc}. An ensemble average of a quantity
can be related to the time average measured in the experiments.

\subsection{Correlation Function\label{In:ssec:corr}}
Thermodynamic properties can be computed by equillibrium statistical
mechanics. On the other hand, {\it Time correlation function} is a
powerful method to understand the evolution of a dynamic system in
non-equillibrium statistical mechanics. Imagine a 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$, at a later time $t_0 + \tau$ this value is
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. Howere, when $\tau$ is large,
this correlation is lost. The correlation function is a measurement of
this relationship and is defined by~\cite{Berne90}
\begin{equation}
C(t) = \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 {\it autocorrelation function}. The defination
of the correlation function for two different variables is similar to
that of autocorrelation function, which is
\begin{equation}
C(t) = \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}.

In section~\ref{In:ssec:average} we know from Eq.~\ref{Ineq:ergodicity}
the relationship between time average and ensemble average. We can put
the correlation function in a classical mechanics form,
\begin{equation}
C(t) = \langle A(0)A(\tau) \rangle = \int d \vec q~^N \int d \vec p~^N A(t) A(t + \tau) \rho(q, p)
\label{Ineq:autocorrelationFunctionCM}
\end{equation}
and
\begin{equation}
C(t) = \langle A(0)B(\tau) \rangle = \int d \vec q~^N \int d \vec p~^N A(t) B(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:	3365
Committed:	Mon Mar 10 21:52:50 2008 UTC (16 years, 6 months ago) by xsun
Content type:	application/x-tex
File size:	30144 byte(s)
Log Message:	writing.