--- trunk/matt_papers/canidacy_paper/canidacy_paper.tex	2002/08/21 21:25:08	95
+++ trunk/matt_papers/canidacy_paper/canidacy_paper.tex	2003/04/28 16:03:10	507
@@ -1,8 +1,12 @@
 \documentclass[11pt]{article}
 
 \usepackage{graphicx}
+\usepackage{color}
+\usepackage{floatflt}
 \usepackage{amsmath}
 \usepackage{amssymb}
+\usepackage{subfigure}
+\usepackage{palatino}
 \usepackage[ref]{overcite}
 
 
@@ -19,6 +23,7 @@
 
 \begin{document}
 
+
 \title{A Mesoscale Model for Phospholipid Simulations}
 
 \author{Matthew A. Meineke\\
@@ -31,63 +36,509 @@ Notre Dame, Indiana 46556}
 
 \section{Background and Research Goals}
 
+Simulations of phospholipid bilayers are, by necessity, quite
+complex. The lipid molecules are large molecules containing many
+atoms, and the head group of the lipid will typically contain charge
+separated ions which set up a large dipole within the molecule. Adding
+to the complexity are the number of water molecules needed to properly
+solvate the lipid bilayer, typically 25 water molecules for every
+lipid molecule. Because of these factors, many current simulations are
+limited in both length and time scale due to to the sheer number of
+calculations performed at every time step and the lifetime of the
+researcher. A typical
+simulation\cite{saiz02,lindahl00,venable00,Marrink01} will have around
+64 phospholipids forming a bilayer approximately 40~$\mbox{\AA}$ by
+50~$\mbox{\AA}$ with roughly 25 waters for every lipid. This means
+there are on the order of 8,000 atoms needed to simulate these systems
+and the trajectories are integrated for times up to 10 ns. 
+
+These limitations make it difficult to study certain biologically
+interesting phenomena that don't fit within the short time and length
+scale requirements. One such phenomenon is the existence of the ripple
+phase ($P_{\beta'}$) of the bilayer between the gel phase
+($L_{\beta'}$) and the fluid phase ($L_{\alpha}$). The $P_{\beta'}$
+phase has been shown to have a ripple period of
+100-200~$\mbox{\AA}$.\cite{katsaras00,sengupta00} A simulation of this
+length scale would require approximately 1,300 lipid molecules and
+roughly 25 waters for every lipid to fully solvate the bilayer. With
+the large number of atoms involved in a simulation of this magnitude,
+steps \emph{must} be taken to simplify the system to the point where
+the numbers of atoms becomes reasonable.
+
+Another system of interest would be drug molecule diffusion through
+the membrane. Due to the fluid-like properties of a lipid membrane,
+not all diffusion takes place at membrane channels. It is of interest
+to study certain molecules that may incorporate themselves directly
+into the membrane. These molecules may then have an appreciable
+waiting time (on the order of nanoseconds) within the
+bilayer. Simulation of such a long time scale again requires
+simplification of the system in order to lower the number of
+calculations needed at each time step or to increase the length of
+each time step.
+
+
 \section{Methodology}
 
-\subsection{Length Scale Simplifications}
+\subsection{Length and Time Scale Simplifications}
 
-The length scale simplifications are aimed at increaseing the number
-of molecules simulated without drastically increasing the
-computational cost of the system. This is done by a combination of
-substituting less expensive interactions for expensive ones and
-decreasing the number of interaction sites per molecule. Namely,
-charge distributions are replaced with dipoles, and unified atoms are
-used in place of water and phospholipid head groups.
+The length scale simplifications are aimed at increasing the number of
+molecules that can be simulated without drastically increasing the
+computational cost of the simulation. This is done through a
+combination of substituting less expensive interactions for expensive
+ones and decreasing the number of interaction sites per
+molecule. Namely, groups of point charges are replaced with single
+point-dipoles, and unified atoms are used in place of water,
+phospholipid head groups, and alkyl groups.
 
-The replacement of charge distributions with dipoles allows us to
-replace an interaction that has a relatively long range, $\frac{1}{r}$
-for the charge charge potential, with that of a relitively short
-range, $\frac{1}{r^{3}}$ for dipole - dipole potentials
-(Equations~\ref{eq:dipolePot} and \ref{eq:chargePot}). This allows us
-to use computaional simplifications algorithms such as Verlet neighbor
-lists,\cite{allen87:csl} which gives computaional scaling by $N$. This
-is in comparison to the Ewald sum\cite{leach01:mm} needed to compute
-the charge - charge interactions which scales at best by $N
-\ln N$.
+The replacement of charges with dipoles allows us to replace an
+interaction that has a relatively long range ($\frac{1}{r}$ for the
+coulomb potential) with that of a relatively short range
+($\frac{1}{r^{3}}$ for dipole - dipole potentials). Combined with
+Verlet neighbor lists,\cite{allen87:csl} this should result in an
+algorithm which scales linearly with increasing system size. This is
+in comparison to the Ewald sum\cite{leach01:mm} needed to compute
+periodic replicas of the coulomb interactions, which scales at
+best\cite{darden93:pme} by $N \ln N$.
 
+The second step taken to simplify the number of calculations is to
+incorporate unified models for groups of atoms. In the case of water,
+we use the soft sticky dipole (SSD) model developed by
+Ichiye\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md}
+(Section~\ref{sec:ssdModel}). For the phospholipids, a unified head
+atom with a dipole will replace the atoms in the head group, while
+unified $\text{CH}_2$ and $\text{CH}_3$ atoms will replace the alkyl
+units in the tails (Section~\ref{sec:lipidModel}). This model is
+similar in practice to that of Lipowsky and Goetz\cite{goetz98}, where
+the whole system is reduced to attractive and repulsive Lennard-Jones
+spheres. However, our model gives a greater level of detail to each
+unified molecule, namely through dipole, bend, and torsion
+interactions.
+
+The time scale simplifications are introduced so that we can take
+longer time steps. By increasing the size of the time steps taken by
+the simulation, we are able to integrate a given length of time using
+fewer calculations. However, care must be taken that any
+simplifications used still conserve the total energy of the
+simulation. In practice, this means taking steps small enough to
+resolve all motion in the system without accidently moving an object
+too far along a repulsive energy surface before it feels the effect of
+the surface.
+
+In our simulation we have chosen to constrain all bonds to be of fixed
+length. This means the bonds are no longer allowed to vibrate about
+their equilibrium positions. Bond vibrations are typically the fastest
+periodic motion in a dynamics simulation. By taking this action, we
+are able to take time steps of 3 fs while still maintaining constant
+energy. This is in contrast to the 1 fs time step typically needed to
+conserve energy when bonds lengths are allowed to oscillate.
+
+\subsection{The Soft Sticky Dipole Water Model}
+\label{sec:ssdModel}
+
+\begin{figure}
+\centering
+\includegraphics[width=50mm]{ssd.epsi}
+\caption{The SSD model with the oxygen and hydrogen atoms drawn in for reference.}
+\label{fig:ssdModel}
+\end{figure} 
+
+The water model used in our simulations is a modified soft
+Stockmayer-sphere model.\cite{stevens95} Like the Stockmayer-sphere, the SSD
+model consists of a Lennard-Jones interaction site and a
+dipole both located at the water's center of mass (Figure
+\ref{fig:ssdModel}). However, the SSD model extends this by adding a
+tetrahedral potential to correct for hydrogen bonding.
+
+The SSD water potential for a pair of water molecules is then given by
+the following equation:
 \begin{equation}
-V^{\text{dp}}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
+V_{\text{SSD}} = V_{\text{LJ}}(r_{i\!j}) + V_{\text{dp}}(\mathbf{r}_{i\!j},
+	\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})
+        + V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
+	\boldsymbol{\Omega}_{j})
+\label{eq:ssdTotPot}
+\end{equation}
+where $\mathbf{r}_{ij}$ is the vector between molecules $i$ and $j$,
+and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ are the
+Euler angles of molecule $i$ or $j$ respectively. $V_{\text{LJ}}$ is
+the Lennard-Jones potential:
+\begin{equation}
+V_{\text{LJ}} = 
+	4\epsilon_{ij} \biggl[
+	\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12}
+	- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6}
+	\biggr]
+\label{eq:lennardJonesPot}
+\end{equation}
+where $\sigma_{ij}$ scales the length of the interaction, and
+$\epsilon_{ij}$ scales the energy of the potential. For SSD,
+$\sigma_{\text{SSD}} = 3.051 \mbox{
+\AA}$ and $\epsilon_{\text{SSD}} = 0.152\text{ kcal/mol}$.
+$V_{\text{dp}}$ is the dipole potential:
+\begin{equation}
+V_{\text{dp}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
 	\boldsymbol{\Omega}_{j}) = \frac{1}{4\pi\epsilon_{0}} \biggl[
 	\frac{\boldsymbol{\mu}_{i} \cdot \boldsymbol{\mu}_{j}}{r^{3}_{ij}}
 	-
-	\frac{3(\boldsymbol{\mu} \cdot \mathbf{r}_{ij}) %
-		(\boldsymbol{\mu} \cdot \mathbf{r}_{ij}) }{r^{5}_{ij}} \biggr]
+	\frac{3(\boldsymbol{\mu}_i \cdot \mathbf{r}_{ij}) %
+		(\boldsymbol{\mu}_j \cdot \mathbf{r}_{ij}) }
+		{r^{5}_{ij}} \biggr]
 \label{eq:dipolePot}
 \end{equation}
+where $\boldsymbol{\mu}_i$ is the dipole vector of molecule $i$,
+$\boldsymbol{\mu}_i$ points along the z-axis in the body fixed
+frame. This frame is then oriented in space by the Euler angles,
+$\boldsymbol{\Omega}_i$. The SSD model specifies a dipole magnitude of
+2.35~D for water.
 
+The hydrogen bonding is modeled by the $V_{\text{sp}}$
+term of the potential. Its form is as follows:
 \begin{equation}
-V^{\text{ch}}_{ij}(\mathbf{r}_{ij}) = \frac{q_{i}q_{j}}%
-	{4\pi\epsilon_{0} r_{ij}}
-\label{eq:chargePot}
-\end{equation} 
+V_{\text{sp}}(\mathbf{r}_{i\!j},\boldsymbol{\Omega}_{i},
+	\boldsymbol{\Omega}_{j}) =
+	v^{\circ}[s(r_{ij})w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
+		\boldsymbol{\Omega}_{j})
+	+
+	s'(r_{ij})w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},
+		\boldsymbol{\Omega}_{j})]
+\label{eq:spPot}
+\end{equation}
+Where $v^\circ$ scales strength of the interaction.
+$w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
+and
+$w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
+are responsible for the tetrahedral potential and a correction to the
+tetrahedral potential respectively. They are,
+\begin{equation}
+w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
+	\sin\theta_{ij} \sin 2\theta_{ij} \cos 2\phi_{ij}
+	+ \sin \theta_{ji} \sin 2\theta_{ji} \cos 2\phi_{ji}
+\label{eq:spPot2}
+\end{equation}
+and
+\begin{equation}
+\begin{split}
+w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) =
+	&(\cos\theta_{ij}-0.6)^2(\cos\theta_{ij} + 0.8)^2 \\
+	&+ (\cos\theta_{ji}-0.6)^2(\cos\theta_{ji} + 0.8)^2 - 2w^{\circ}
+\end{split}
+\label{eq:spCorrection}
+\end{equation}
+The angles $\theta_{ij}$ and $\phi_{ij}$ are defined by the spherical
+coordinates of the position of molecule $j$ in the reference frame
+fixed on molecule $i$ with the z-axis aligned with the dipole moment.
+The correction
+$w^{x}_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
+is needed because
+$w_{ij}(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j})$
+vanishes when $\theta_{ij}$ is $0^\circ$ or $180^\circ$.
 
-The second step taken to simplify the number of calculationsis to
-incorporate unified models for groups of atoms. In the case of water,
-we use the soft sticky dipole (SSD) model developed by
-Ichiye\cite{Liu96} (Section~\ref{sec:ssdModel}). For the phospholipids, a
-unified head atom with a dipole will replace the atoms in the head
-group, while unified $\text{CH}_2$ and $\text{CH}_3$ atoms will
-replace the alkanes in the tails (Section~\ref{sec:lipidModel}).
+Finally, the sticky potential is scaled by a switching function,
+$s(r_{ij})$, that scales smoothly between 0 and 1. It is represented
+by:
+\begin{equation}
+s(r_{ij}) = 
+	\begin{cases}
+	1&	\text{if $r_{ij} < r_{L}$}, \\
+	\frac{(r_{U} - r_{ij})^2 (r_{U} + 2r_{ij} 
+		- 3r_{L})}{(r_{U}-r_{L})^3}&
+		\text{if $r_{L} \leq r_{ij} \leq r_{U}$},\\
+	0&	\text{if $r_{ij} \geq r_{U}$}.
+	\end{cases}
+\label{eq:spCutoff}
+\end{equation}
 
+Despite the apparent complexity of Equation \ref{eq:spPot}, the SSD
+model is still computationally inexpensive. This is due to Equation
+\ref{eq:spCutoff}. With $r_{L}$ being 2.75~$\mbox{\AA}$ and $r_{U}$
+being equal to either 3.35~$\mbox{\AA}$ for $s(r_{ij})$ or
+4.0~$\mbox{\AA}$ for $s'(r_{ij})$, the sticky potential is only active
+over an extremely short range, and then only with other SSD
+molecules. Therefore, it's predominant interaction is through the
+point dipole and the Lennard-Jones sphere. Of these, only the dipole
+interaction can be considered ``long-range''.
 
-\subsection{Time Scale Simplifications}
-
-\subsection{The Soft Sticky Water Model}
-\label{sec:ssdModel}
-
 \subsection{The Phospholipid Model}
 \label{sec:lipidModel}
 
+\begin{figure}
+\centering
+\includegraphics[angle=-90,width=80mm]{lipidModel.epsi}
+\caption{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ is the bend angle, $\mu$ is the dipole moment of the head group, and n is the chain length.}
+\label{fig:lipidModel}
+\end{figure}
 
-\bibliographystyle{achemso}
-\bibliography{canidacy_paper} 
-\end{document}
\ No newline at end of file
+The lipid molecules in our simulations are unified atom models. Figure
+\ref{fig:lipidModel} shows a schematic for one of our
+lipids. The head group of the phospholipid is replaced by a single
+Lennard-Jones sphere with a freely oriented dipole placed at it's
+center. The magnitude of the dipole moment is 20.6 D, chosen to match
+that of DPPC\cite{Cevc87}. The tail atoms are unified $\text{CH}_2$
+and $\text{CH}_3$ atoms and are also modeled as Lennard-Jones
+spheres. The total potential for the lipid is represented by Equation
+\ref{eq:lipidModelPot}.
+
+\begin{equation}
+V_{\text{lipid}} =
+	\sum_{i}V_{i}^{\text{internal}}
+	+ \sum_i \sum_{j>i} \sum_{\alpha_i} 
+	\sum_{\beta_j}
+	V_{\text{LJ}}(r_{\alpha_{i}\beta_{j}})
+	+\sum_i\sum_{j>i}V_{\text{dp}}(r_{1_i,1_j},\Omega_{1_i},\Omega_{1_j})
+\label{eq:lipidModelPot}
+\end{equation}
+where,
+\begin{equation}
+V_{i}^{\text{internal}} = 
+	\sum_{\text{bends}}V_{\text{bend}}(\theta_{\alpha\beta\gamma})
+	+ \sum_{\text{torsions}}V_{\text{tors.}}(\phi_{\alpha\beta\gamma\zeta})
+	+ \sum_{\alpha_i} \sum_{\beta_i > (\alpha_i + 4)}V_{\text{LJ}}
+		(r_{\alpha_i \beta_i})
+\label{eq:lipidModelPotInternal}
+\end{equation}
+
+The non-bonded interactions, $V_{\text{LJ}}$ and $V_{\text{dp}}$, are
+the Lennard-Jones and dipole-dipole interactions respectively. For the
+bonded potentials, only the bend and the torsional potentials are
+calculated. The bond potential is not calculated, and the bond lengths
+are constrained via RATTLE.\cite{leach01:mm} The bend potential is of
+the form: 
+\begin{equation}
+V_{\text{bend}}(\theta_{\alpha\beta\gamma}) 
+	= k_{\theta}\frac{(\theta_{\alpha\beta\gamma} - \theta_0)^2}{2}
+\label{eq:bendPot}
+\end{equation}
+Where $k_{\theta}$ sets the stiffness of the bend potential, and $\theta_0$
+sets the equilibrium bend angle. The torsion potential was given by:
+\begin{equation}
+V_{\text{tors.}}(\phi_{\alpha\beta\gamma\zeta})
+	= c_1 [1+\cos\phi_{\alpha\beta\gamma\zeta}]
+	+ c_2 [1 - \cos(2\phi_{\alpha\beta\gamma\zeta})] 
+	+ c_3 [1 + \cos(3\phi_{\alpha\beta\gamma\zeta})]
+\label{eq:torsPot}
+\end{equation}
+All parameters for bonded and non-bonded potentials in the tail atoms
+were taken from TraPPE.\cite{Siepmann1998} The bonded interactions for
+the head atom were also taken from TraPPE, however it's dipole moment
+and mass were based on the properties of the phosphatidylcholine head
+group. The Lennard-Jones parameter for the head group was chosen such
+that it was roughly twice the size of a $\text{CH}_3$ atom, and it's
+well depth was set to be approximately equal to that of $\text{CH}_3$.
+
+\section{Initial Simulation: 25 lipids in water}
+\label{sec:5x5}
+
+\subsection{Starting Configuration and Parameters}
+\label{sec:5x5Start}
+
+\begin{figure}
+\centering
+\mbox{
+	\subfigure[The starting configuration of the 25 lipid system.]{%
+		\label{fig:5x5Start}%
+		\includegraphics[width=70mm]{5x5-initial.eps}}\quad
+	\subfigure[The 25 lipid system at 6.27~ns.]{%
+		\label{fig:5x5Final}%
+		\includegraphics[width=70mm]{5x5-6.27ns.eps}}
+}
+\caption{Snapshots of the 25 lipid system. Carbon tail atoms are drawn in gray, the phospholipid head groups are colored blue, and the waters are scaled down for visibility. A box has been drawn around the periodic image.}
+\end{figure}
+
+Our first simulation is an array of 25 single chain lipids in a sea
+of water (Figure \ref{fig:5x5Start}). The total number of water
+molecules is 1386, giving a final of water concentration of 70\%
+wt. The simulation box measures 34.5~$\mbox{\AA}$ x 39.4~$\mbox{\AA}$
+x 39.4~$\mbox{\AA}$ with periodic boundary conditions imposed. The
+system is simulated in the micro-canonical (NVE) ensemble with an
+average temperature of 300~K.
+
+\subsection{Results}
+\label{sec:5x5Results}
+
+\begin{figure}
+\centering
+	\subfigure[The self correlation of the phospholipid head groups. $g(r)$ is on the top, the bottom chart is the $g_\gamma(r)$.]{%
+		\label{fig:5x5HHCorr}%
+		\includegraphics[angle=-90,width=80mm]{all5x5-HEAD-HEAD.epsi}%
+		}
+	\subfigure[The $g(r)$ for the $\text{CH}_2$ molecules in the chain tails]{%
+		\label{fig:5x5CCg}%
+		\includegraphics[angle=-90,width=80mm]{all5x5-CH2-CH2.epsi}}
+	\subfigure%
+	[The pair correlations between the head groups and the water]{%
+		\label{fig:5x5HXCorr}%
+		\includegraphics[angle=-90,width=80mm]{all5x5-HEAD-X.epsi}}
+\caption{The pair correlation functions for the 25 lipid system}
+\end{figure}
+		
+
+Figure \ref{fig:5x5Final} shows a snapshot of the system at
+6.27~ns. Note that the system has spontaneously self assembled into a
+bilayer. Discussion of the length scales of the bilayer will follow in
+this section. However, it is interesting to note a key qualitative
+property of the system revealed by this snapshot, the tail chains are
+tilted to the bilayer normal. This is usually indicative of the gel
+($L_{\beta'}$) phase. In this system, the box size is probably too
+small for the bilayer to relax to the fluid ($P_{\alpha}$) phase. This
+demonstrates a need for an isobaric-isothermal ensemble where the box
+size may relax or expand to keep the system at 1~atm.
+
+The simulation was analyzed using the radial distribution function,
+$g(r)$, which has the form:
+\begin{equation}
+g(r) = \frac{V}{N_{\text{pairs}}}\langle \sum_{i} \sum_{j > i} 
+	\delta(|\mathbf{r} - \mathbf{r}_{ij}|) \rangle
+\label{eq:gofr}
+\end{equation}
+Equation \ref{eq:gofr} gives us information about the spacing of two
+species as a function of radius. Essentially, if the observer were
+located at atom $i$ and were looking out in all directions, $g(r)$
+shows the relative density of atom $j$ at any given radius, $r$,
+normalized by the expected density of atom $j$ at $r$. In a
+homogeneously mixed fluid, $g(r)$ will approach 1 at large $r$, as a
+fluid contains no long range structure to contribute peaks in the
+number density.
+
+For the species containing dipoles, a second pair-wise distribution
+function was used, $g_{\gamma}(r)$. It is of the form:
+\begin{equation}
+g_{\gamma}(r) = \langle \sum_i \sum_{j>i}
+	(\cos \gamma_{ij}) \delta(| \mathbf{r} - \mathbf{r}_{ij}|) \rangle
+\label{eq:gammaofr}
+\end{equation}
+Where $\gamma_{ij}$ is the angle between the dipole of atom $j$ with
+respect to the dipole of atom $i$. This correlation will vary between
+$+1$ and $-1$ when the two dipoles are perfectly aligned and
+anti-aligned respectively. This then gives us information about how
+directional species are aligned with each other as a function of
+distance.
+
+Figure \ref{fig:5x5HHCorr} shows the two self correlation functions
+for the Head groups of the lipids. The first peak in the $g(r)$ at
+4.03~$\mbox{\AA}$ is the nearest neighbor separation of the heads of
+two lipids. This corresponds to a maximum in the $g_{\gamma}(r)$ which
+means that the two neighbors on the same leaf have their dipoles
+aligned. The broad peak at 6.5~$\mbox{\AA}$ is the inter-surface
+spacing. Here, there is a corresponding anti-alignment in the angular
+correlation. This means that although the dipoles are aligned on the
+same monolayer, the dipoles will orient themselves to be anti-aligned
+on a opposite facing monolayer. With this information, the two peaks
+at 7.0~$\mbox{\AA}$ and 7.4~$\mbox{\AA}$ are head groups on the same
+monolayer, and they are the second nearest neighbors to the head
+group. The peak is likely a split peak because of the small statistics
+of this system. Finally, the peak at 8.0~$\mbox{\AA}$ is likely the
+second nearest neighbor on the opposite monolayer because of the
+anti-alignment evident in the angular correlation.
+
+Figure \ref{fig:5x5CCg} shows the radial distribution function for the
+$\text{CH}_2$ unified atoms. The spacing of the atoms along the tail
+chains accounts for the regularly spaced sharp peaks, but the broad
+underlying peak with its maximum at 4.6~$\mbox{\AA}$ is the
+distribution of chain-chain distances between two lipids. The final
+Figure, Figure \ref{fig:5x5HXCorr}, includes the correlation functions
+between the Head group and the SSD atoms. The peak in $g(r)$ at
+3.6~$\mbox{\AA}$ is the distance of closest approach between the two,
+and $g_{\gamma}(r)$ shows that the SSD atoms will align their dipoles
+with the head groups at close distance. However, as one increases the
+distance, the SSD atoms are no longer aligned.
+
+\section{Second Simulation: 50 randomly oriented lipids in water}
+\label{sec:r50}
+
+\subsection{Starting Configuration and Parameters}
+\label{sec:r50Start}
+
+\begin{figure}
+\centering
+\mbox{
+	\subfigure[The starting configuration of the 50 lipid system.]{%
+		\label{fig:r50Start}%
+		\includegraphics[width=70mm]{r50-initial.eps}}\quad	
+	\subfigure[The 50 lipid system at 2.21~ns]{%
+		\label{fig:r50Final}%	
+		\includegraphics[width=70mm]{r50-2.21ns.eps}}
+}
+\caption{Snapshots of the 50 lipid system}
+\end{figure}
+
+The second simulation consists of 50 single chained lipid molecules
+embedded in a sea of 1384 SSD waters (54\% wt.). The lipids in this
+simulation were started with random orientation and location (Figure
+\ref{fig:r50Start} ) The simulation box measured 26.6~$\mbox{\AA}$ x
+26.6~$\mbox{\AA}$ x 108.4~$\mbox{\AA}$ with periodic boundary conditions
+imposed. The simulation was run in the NVE ensemble with an average
+temperature of 300~K.
+
+\subsection{Results}
+\label{sec:r50Results}
+
+\begin{figure}
+\centering
+	\subfigure[The self correlation of the phospholipid head groups.]{%
+		\label{fig:r50HHCorr}%
+		\includegraphics[angle=-90,width=80mm]{r50-HEAD-HEAD.epsi}%
+		}
+	\subfigure%
+	[The pair correlations between the head groups and the water]{%
+		\label{fig:r50HXCorr}%
+		\includegraphics[angle=-90,width=80mm]{r50-HEAD-X.epsi}%
+		}
+	\subfigure[The $g(r)$ for the $\text{CH}_2$ molecules in the chain tails]{%
+		\label{fig:r50CCg}%
+		\includegraphics[angle=-90,width=80mm]{r50-CH2-CH2.epsi}}
+
+\caption{The pair correlation functions for the 50 lipid system}
+\end{figure}
+
+Figure \ref{fig:r50Final} is a snapshot of the system at 2.21~ns. Here
+we see that the system has already aggregated into several micelles
+and two are even starting to merge. It will be interesting to watch as
+this simulation continues what the total time scale for the micelle
+aggregation and bilayer formation will be, in Marrink's\cite{Marrink01}
+simulation, bilayer aggregation is predicted to occur around 10~ns.
+
+Figures \ref{fig:r50HHCorr}, \ref{fig:r50HXCorr}, and \ref{fig:r50CCg} are
+the same correlation functions for the random 50 simulation as for the
+previous simulation of 25 lipids. What is most interesting to note, is
+the high degree of similarity between the correlation functions
+between systems. Even though the 25 lipid simulation formed a bilayer
+and the random 50 simulation is still in the micelle stage, both have
+an inter-surface spacing of 6.5~$\mbox{\AA}$ with the same
+characteristic anti-alignment between surfaces. Not as surprising, is
+the consistency of the closest packing statistics between
+systems. Namely, a head-head closest approach distance of
+4~$\mbox{\AA}$, and similar findings for the chain-chain and
+head-water distributions as in the 25 lipid system.
+
+\section{Future Directions}
+
+Current simulations indicate that our model is a feasible one, however
+improvements will need to be made to allow the system to be simulated
+in the isobaric-isothermal ensemble. This will relax the system to an
+equilibrium configuration at room temperature and pressure allowing us
+to compare our model to experimental results. Also, we are in the
+process of parallelizing the code for an even greater speedup. This
+will allow us to simulate large enough systems to examine phenomena
+such as the ripple phase and drug molecule diffusion
+
+Once the work has been completed on the simulation engine, we will
+then use it to explore the phase diagram for our model. By
+characterizing how our model parameters affect the bilayer properties,
+we will tailor our model to more closely match real biological
+molecules. With this information, we will then incorporate
+biologically relevant molecules into the system and observe their
+transport properties across the membrane.
+
+\section{Acknowledgments}
+
+I would like to thank Dr.~J.~Daniel Gezelter for his guidance on this
+project. I would also like to acknowledge the following for their help
+and discussions during this project: Christopher Fennell, Charles
+Vardeman, Teng Lin, Megan Sprague, Patrick Conforti, and Dan
+Combest. Funding for this project came from the National Science
+Foundation.
+
+\pagebreak 
+\bibliographystyle{achemso} 
+\bibliography{canidacy_paper}
+\end{document}