--- trunk/mattDisertation/Introduction.tex	2004/01/18 05:26:10	955
+++ trunk/mattDisertation/Introduction.tex	2004/01/22 21:13:55	977
@@ -53,10 +53,10 @@ The equation can be recast as:
 \end{equation}
 The equation can be recast as:
 \begin{equation}
-I = (b-a)<f(x)>
+I = (b-a)\langle f(x) \rangle
 \label{eq:MCex2}
 \end{equation}
-Where $<f(x)>$ is the unweighted average over the interval
+Where $\langle f(x) \rangle$ is the unweighted average over the interval
 $[a,b]$. The calculation of the integral could then be solved by
 randomly choosing points along the interval $[a,b]$ and calculating
 the value of $f(x)$ at each point. The accumulated average would then
@@ -66,16 +66,18 @@ integrals of the form:
 However, in Statistical Mechanics, one is typically interested in
 integrals of the form:
 \begin{equation}
-<A> = \frac{A}{exp^{-\beta}}
+\langle A \rangle = \frac{\int d^N \mathbf{r}~A(\mathbf{r}^N)% 
+	e^{-\beta V(\mathbf{r}^N)}}%
+	{\int d^N \mathbf{r}~e^{-\beta V(\mathbf{r}^N)}}
 \label{eq:mcEnsAvg}
 \end{equation}
-Where $r^N$ stands for the coordinates of all $N$ particles and $A$ is
-some observable that is only dependent on position. $<A>$ is the
-ensemble average of $A$ as presented in
-Sec.~\ref{introSec:statThermo}. Because $A$ is independent of
-momentum, the momenta contribution of the integral can be factored
-out, leaving the configurational integral. Application of the brute
-force method to this system would yield highly inefficient
+Where $\mathbf{r}^N$ stands for the coordinates of all $N$ particles
+and $A$ is some observable that is only dependent on
+position. $\langle A \rangle$ is the ensemble average of $A$ as
+presented in Sec.~\ref{introSec:statThermo}. Because $A$ is
+independent of momentum, the momenta contribution of the integral can
+be factored out, leaving the configurational integral. Application of
+the brute force method to this system would yield highly inefficient
 results. Due to the Boltzman weighting of this integral, most random
 configurations will have a near zero contribution to the ensemble
 average. This is where a importance sampling comes into
@@ -85,23 +87,258 @@ Eq.~\ref{eq:MCex1} rewritten to be:
 which the random configurations are chosen in order to more
 efficiently calculate the integral.\cite{Frenkel1996} Consider again
 Eq.~\ref{eq:MCex1} rewritten to be:
+\begin{equation}
+I = \int^b_a \frac{f(x)}{\rho(x)} \rho(x) dx
+\label{introEq:Importance1}
+\end{equation}
+Where $\rho(x)$ is an arbitrary probability distribution in $x$.  If
+one conducts $\tau$ trials selecting a random number, $\zeta_\tau$,
+from the distribution $\rho(x)$ on the interval $[a,b]$, then
+Eq.~\ref{introEq:Importance1} becomes
+\begin{equation}
+I= \biggl \langle \frac{f(x)}{\rho(x)} \biggr \rangle_{\text{trials}}
+\label{introEq:Importance2}
+\end{equation}
+Looking at Eq.~\ref{eq:mcEnsAvg}, and realizing
+\begin {equation}
+\rho_{kT}(\mathbf{r}^N) = 
+	\frac{e^{-\beta V(\mathbf{r}^N)}}
+	{\int d^N \mathbf{r}~e^{-\beta V(\mathbf{r}^N)}}
+\label{introEq:MCboltzman}
+\end{equation}
+Where $\rho_{kT}$ is the boltzman distribution.  The ensemble average
+can be rewritten as
+\begin{equation}
+\langle A \rangle = \int d^N \mathbf{r}~A(\mathbf{r}^N) 
+	\rho_{kT}(\mathbf{r}^N)
+\label{introEq:Importance3}
+\end{equation}
+Applying Eq.~\ref{introEq:Importance1} one obtains
+\begin{equation}
+\langle A \rangle = \biggl \langle
+	\frac{ A \rho_{kT}(\mathbf{r}^N) }
+	{\rho(\mathbf{r}^N)} \biggr \rangle_{\text{trials}}
+\label{introEq:Importance4}
+\end{equation}
+By selecting $\rho(\mathbf{r}^N)$ to be $\rho_{kT}(\mathbf{r}^N)$
+Eq.~\ref{introEq:Importance4} becomes
+\begin{equation}
+\langle A \rangle = \langle A(\mathbf{r}^N) \rangle_{\text{trials}}
+\label{introEq:Importance5}
+\end{equation}
+The difficulty is selecting points $\mathbf{r}^N$ such that they are
+sampled from the distribution $\rho_{kT}(\mathbf{r}^N)$.  A solution
+was proposed by Metropolis et al.\cite{metropolis:1953} which involved
+the use of a Markov chain whose limiting distribution was
+$\rho_{kT}(\mathbf{r}^N)$.
 
+\subsubsection{\label{introSec:markovChains}Markov Chains}
 
+A Markov chain is a chain of states satisfying the following
+conditions:\cite{leach01:mm}
+\begin{enumerate}
+\item The outcome of each trial depends only on the outcome of the previous trial.
+\item Each trial belongs to a finite set of outcomes called the state space.
+\end{enumerate}
+If given two configuartions, $\mathbf{r}^N_m$ and $\mathbf{r}^N_n$,
+$\rho_m$ and $\rho_n$ are the probablilities of being in state
+$\mathbf{r}^N_m$ and $\mathbf{r}^N_n$ respectively.  Further, the two
+states are linked by a transition probability, $\pi_{mn}$, which is the
+probability of going from state $m$ to state $n$.
 
-\subsection{\label{introSec:md}Molecular Dynamics Simulations}
+\newcommand{\accMe}{\operatorname{acc}}
 
-time averages
+The transition probability is given by the following:
+\begin{equation}
+\pi_{mn} = \alpha_{mn} \times \accMe(m \rightarrow n)
+\label{introEq:MCpi}
+\end{equation}
+Where $\alpha_{mn}$ is the probability of attempting the move $m
+\rightarrow n$, and $\accMe$ is the probability of accepting the move
+$m \rightarrow n$.  Defining a probability vector,
+$\boldsymbol{\rho}$, such that
+\begin{equation}
+\boldsymbol{\rho} = \{\rho_1, \rho_2, \ldots \rho_m, \rho_n, 
+	\ldots \rho_N \}
+\label{introEq:MCrhoVector}
+\end{equation}
+a transition matrix $\boldsymbol{\Pi}$ can be defined,
+whose elements are $\pi_{mn}$, for each given transition.  The
+limiting distribution of the Markov chain can then be found by
+applying the transition matrix an infinite number of times to the
+distribution vector.
+\begin{equation}
+\boldsymbol{\rho}_{\text{limit}} = 
+	\lim_{N \rightarrow \infty} \boldsymbol{\rho}_{\text{initial}}
+	\boldsymbol{\Pi}^N
+\label{introEq:MCmarkovLimit}
+\end{equation}
+The limiting distribution of the chain is independent of the starting
+distribution, and successive applications of the transition matrix
+will only yield the limiting distribution again.
+\begin{equation}
+\boldsymbol{\rho}_{\text{limit}} = \boldsymbol{\rho}_{\text{initial}}
+	\boldsymbol{\Pi}
+\label{introEq:MCmarkovEquil}
+\end{equation}
 
-time integrating schemes
+\subsubsection{\label{introSec:metropolisMethod}The Metropolis Method}
 
-time reversible
+In the Metropolis method\cite{metropolis:1953}
+Eq.~\ref{introEq:MCmarkovEquil} is solved such that
+$\boldsymbol{\rho}_{\text{limit}}$ matches the Boltzman distribution
+of states.  The method accomplishes this by imposing the strong
+condition of microscopic reversibility on the equilibrium
+distribution.  Meaning, that at equilibrium the probability of going
+from $m$ to $n$ is the same as going from $n$ to $m$.
+\begin{equation}
+\rho_m\pi_{mn} = \rho_n\pi_{nm}
+\label{introEq:MCmicroReverse}
+\end{equation}
+Further, $\boldsymbol{\alpha}$ is chosen to be a symetric matrix in
+the Metropolis method.  Using Eq.~\ref{introEq:MCpi},
+Eq.~\ref{introEq:MCmicroReverse} becomes
+\begin{equation}
+\frac{\accMe(m \rightarrow n)}{\accMe(n \rightarrow m)} =
+	\frac{\rho_n}{\rho_m}
+\label{introEq:MCmicro2}
+\end{equation}
+For a Boltxman limiting distribution,
+\begin{equation}
+\frac{\rho_n}{\rho_m} = e^{-\beta[\mathcal{U}(n) - \mathcal{U}(m)]}
+	= e^{-\beta \Delta \mathcal{U}}
+\label{introEq:MCmicro3}
+\end{equation}
+This allows for the following set of acceptance rules be defined:
+\begin{equation}
+EQ Here
+\end{equation}
 
-symplectic methods
+Using the acceptance criteria from Eq.~\ref{fix} the Metropolis method
+proceeds as follows
+\begin{itemize}
+\item Generate an initial configuration $fix$ which has some finite probability in $fix$.
+\item Modify $fix$, to generate configuratioon $fix$.
+\item If configuration $n$ lowers the energy of the system, accept the move with unity ($fix$ becomes $fix$).  Otherwise accept with probability $fix$.
+\item Accumulate the average for the configurational observable of intereest.
+\item Repeat from step 2 until average converges.
+\end{itemize}
+One important note is that the average is accumulated whether the move
+is accepted or not, this ensures proper weighting of the average.
+Using Eq.~\ref{fix} it becomes clear that the accumulated averages are
+the ensemble averages, as this method ensures that the limiting
+distribution is the Boltzman distribution.
 
-Extended ensembles (NVT NPT)
+\subsection{\label{introSec:MD}Molecular Dynamics Simulations}
 
-constrained dynamics
+The main simulation tool used in this research is Molecular Dynamics.
+Molecular Dynamics is when the equations of motion for a system are
+integrated in order to obtain information about both the positions and
+momentum of a system, allowing the calculation of not only
+configurational observables, but momenta dependent ones as well:
+diffusion constants, velocity auto correlations, folding/unfolding
+events, etc.  Due to the principle of ergodicity, Eq.~\ref{fix}, the
+average of these observables over the time period of the simulation
+are taken to be the ensemble averages for the system.
 
+The choice of when to use molecular dynamics over Monte Carlo
+techniques, is normally decided by the observables in which the
+researcher is interested.  If the observabvles depend on momenta in
+any fashion, then the only choice is molecular dynamics in some form.
+However, when the observable is dependent only on the configuration,
+then most of the time Monte Carlo techniques will be more efficent.
+
+The focus of research in the second half of this dissertation is
+centered around the dynamic properties of phospholipid bilayers,
+making molecular dynamics key in the simulation of those properties.
+
+\subsubsection{Molecular dynamics Algorithm}
+
+To illustrate how the molecular dynamics technique is applied, the
+following sections will describe the sequence involved in a
+simulation.  Sec.~\ref{fix} deals with the initialization of a
+simulation.  Sec.~\ref{fix} discusses issues involved with the
+calculation of the forces.  Sec.~\ref{fix} concludes the algorithm
+discussion with the integration of the equations of motion. \cite{fix}
+
+\subsubsection{initialization}
+
+When selecting the initial configuration for the simulation it is
+important to consider what dynamics one is hoping to observe.
+Ch.~\ref{fix} deals with the formation and equilibrium dynamics of
+phospholipid membranes.  Therefore in these simulations initial
+positions were selected that in some cases dispersed the lipids in
+water, and in other cases structured the lipids into preformed
+bilayers.  Important considerations at this stage of the simulation are:
+\begin{itemize}
+\item There are no major overlaps of molecular or atomic orbitals
+\item Velocities are chosen in such a way as to not gie the system a non=zero total momentum or angular momentum.
+\item It is also sometimes desireable to select the velocities to correctly sample the target temperature.
+\end{itemize}
+
+The first point is important due to the amount of potential energy
+generated by having two particles too close together.  If overlap
+occurs, the first evaluation of forces will return numbers so large as
+to render the numerical integration of teh motion meaningless.  The
+second consideration keeps the system from drifting or rotating as a
+whole.  This arises from the fact that most simulations are of systems
+in equilibrium in the absence of outside forces.  Therefore any net
+movement would be unphysical and an artifact of the simulation method
+used.  The final point addresses teh selection of the magnitude of the
+initial velocities.  For many simulations it is convienient to use
+this opportunity to scale the amount of kinetic energy to reflect the
+desired thermal distribution of the system.  However, it must be noted
+that most systems will require further velocity rescaling after the
+first few initial simulation steps due to either loss or gain of
+kinetic energy from energy stored in potential degrees of freedom.
+
+\subsubsection{Force Evaluation}
+
+The evaluation of forces is the most computationally expensive portion
+of a given molecular dynamics simulation.  This is due entirely to the
+evaluation of long range forces in a simulation, typically pair-wise.
+These forces are most commonly the Van der Waals force, and sometimes
+Coulombic forces as well.  For a pair-wise force, there are $fix$
+pairs to be evaluated, where $n$ is the number of particles in the
+system.  This leads to the calculations scaling as $fix$, making large
+simulations prohibitive in the absence of any computation saving
+techniques.
+
+Another consideration one must resolve, is that in a given simulation
+a disproportionate number of the particles will feel the effects of
+the surface. \cite{fix} For a cubic system of 1000 particles arranged
+in a $10x10x10$ cube, 488 particles will be exposed to the surface.
+Unless one is simulating an isolated particle group in a vacuum, the
+behavior of the system will be far from the desired bulk
+charecteristics.  To offset this, simulations employ the use of
+periodic boundary images. \cite{fix}
+
+The technique involves the use of an algorithm that replicates the
+simulation box on an infinite lattice in cartesian space.  Any given
+particle leaving the simulation box on one side will have an image of
+itself enter on the opposite side (see Fig.~\ref{fix}).
+\begin{equation}
+EQ Here
+\end{equation}
+In addition, this sets that any given particle pair has an image, real
+or periodic, within $fix$ of each other.  A discussion of the method
+used to calculate the periodic image can be found in Sec.\ref{fix}.
+
+Returning to the topic of the computational scale of the force
+evaluation, the use of periodic boundary conditions requires that a
+cutoff radius be employed.  Using a cutoff radius improves the
+efficiency of the force evaluation, as particles farther than a
+predetermined distance, $fix$, are not included in the
+calculation. \cite{fix} In a simultation with periodic images, $fix$
+has a maximum value of $fix$.  Fig.~\ref{fix} illustrates how using an
+$fix$ larger than this value, or in the extreme limit of no $fix$ at
+all, the corners of the simulation box are unequally weighted due to
+the lack of particle images in the $x$, $y$, or $z$ directions past a
+disance of $fix$.
+
+With the use of an $fix$, however, comes a discontinuity in the potential energy curve (Fig.~\ref{fix}).
+
+
 \section{\label{introSec:chapterLayout}Chapter Layout}
 
 \subsection{\label{introSec:RSA}Random Sequential Adsorption}