--- trunk/mattDisertation/Introduction.tex	2004/02/03 17:41:56	1008
+++ trunk/mattDisertation/Introduction.tex	2004/02/03 21:30:23	1014
@@ -285,7 +285,7 @@ sampled from the distribution $\rho_{kT}(\mathbf{r}^N)
 \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
+was proposed by Metropolis \emph{et al}.\cite{metropolis:1953} which involved
 the use of a Markov chain whose limiting distribution was
 $\rho_{kT}(\mathbf{r}^N)$.
 
@@ -481,9 +481,9 @@ itself enter on the opposite side (see Fig.~\ref{intro
 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{introFig:pbc}).  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
+addition, this sets that any two particles have an image, real or
+periodic, within $\text{box}/2$ of each other.  A discussion of the
+method used to calculate the periodic image can be found in
 Sec.\ref{oopseSec:pbc}.
 
 \begin{figure}
@@ -557,15 +557,19 @@ q(t-\Delta t)= q(t) - v(t)\Delta t + \frac{F(t)}{2m}\D
 	\mathcal{O}(\Delta t^4) 
 \label{introEq:verletBack}
 \end{equation}
-Adding together Eq.~\ref{introEq:verletForward} and
+Where $m$ is the mass of the particle, $q(t)$ is the position at time
+$t$, $v(t)$ the velocity, and $F(t)$ the force acting on the
+particle. Adding together Eq.~\ref{introEq:verletForward} and
 Eq.~\ref{introEq:verletBack} results in,
 \begin{equation}
-eq here
+q(t+\Delta t)+q(t-\Delta t) = 
+	2q(t) + \frac{F(t)}{m}\Delta t^2 + \mathcal{O}(\Delta t^4) 
 \label{introEq:verletSum}
 \end{equation}
 Or equivalently,
 \begin{equation}
-eq here
+q(t+\Delta t) \approx 
+	2q(t) - q(t-\Delta t) + \frac{F(t)}{m}\Delta t^2
 \label{introEq:verletFinal}
 \end{equation}
 Which contains an error in the estimate of the new positions on the
@@ -573,14 +577,13 @@ with a velocity reformulation of the Verlet method.\ci
 
 In practice, however, the simulations in this research were integrated
 with a velocity reformulation of the Verlet method.\cite{allen87:csl}
-\begin{equation}
-eq here
-\label{introEq:MDvelVerletPos}
-\end{equation}
-\begin{equation}
-eq here
+\begin{align}
+q(t+\Delta t) &= q(t) + v(t)\Delta t + \frac{F(t)}{2m}\Delta t^2 %
+\label{introEq:MDvelVerletPos} \\%
+%
+v(t+\Delta t) &= v(t) + \frac{\Delta t}{2m}[F(t) + F(t+\Delta t)] %
 \label{introEq:MDvelVerletVel}
-\end{equation}
+\end{align}
 The original Verlet algorithm can be regained by substituting the
 velocity back into Eq.~\ref{introEq:MDvelVerletPos}.  The Verlet
 formulations are chosen in this research because the algorithms have
@@ -602,74 +605,103 @@ ensemble average of the observable being measured.  Fr
 reversible.  The fact that it shadows the true Hamiltonian in phase
 space is acceptable in actual simulations as one is interested in the
 ensemble average of the observable being measured.  From the ergodic
-hypothesis (Sec.~\ref{introSec:StatThermo}), it is known that the time
+hypothesis (Sec.~\ref{introSec:statThermo}), it is known that the time
 average will match the ensemble average, therefore two similar
 trajectories in phase space should give matching statistical averages.
 
 \subsection{\label{introSec:MDfurther}Further Considerations}
+
 In the simulations presented in this research, a few additional
 parameters are needed to describe the motions.  The simulations
-involving water and phospholipids in Ch.~\ref{chaptLipids} are
+involving water and phospholipids in Ch.~\ref{chapt:lipid} are
 required to integrate the equations of motions for dipoles on atoms.
 This involves an additional three parameters be specified for each
 dipole atom: $\phi$, $\theta$, and $\psi$.  These three angles are
 taken to be the Euler angles, where $\phi$ is a rotation about the
 $z$-axis, and $\theta$ is a rotation about the new $x$-axis, and
 $\psi$ is a final rotation about the new $z$-axis (see
-Fig.~\ref{introFig:euleerAngles}).  This sequence of rotations can be
-accumulated into a single $3 \times 3$ matrix $\mathbf{A}$
+Fig.~\ref{introFig:eulerAngles}).  This sequence of rotations can be
+accumulated into a single $3 \times 3$ matrix, $\mathbf{A}$,
 defined as follows:
 \begin{equation}
-eq here
+\mathbf{A} = 
+\begin{bmatrix}
+	\cos\phi\cos\psi-\sin\phi\cos\theta\sin\psi &%
+	\sin\phi\cos\psi+\cos\phi\cos\theta\sin\psi &%
+	\sin\theta\sin\psi \\%
+	%
+	-\cos\phi\sin\psi-\sin\phi\cos\theta\cos\psi &%
+	-\sin\phi\sin\psi+\cos\phi\cos\theta\cos\psi &%
+	\sin\theta\cos\psi \\%
+	%
+	\sin\phi\sin\theta &%
+	-\cos\phi\sin\theta &%
+	\cos\theta
+\end{bmatrix}
 \label{introEq:EulerRotMat}
 \end{equation}
 
-The equations of motion for Euler angles can be written down as
-\cite{allen87:csl}
-\begin{equation}
-eq here
-\label{introEq:MDeuleeerPsi}
-\end{equation}
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{eulerRotFig.eps}
+\caption[Euler rotation of Cartesian coordinates]{The rotation scheme for Euler angles. First is a rotation of $\phi$ about the $z$ axis (blue rotation). Next is a rotation of $\theta$ about the new $x\prime$ axis (green rotation). Lastly is a final rotation of $\psi$ about the new $z\prime$ axis (red rotation).}
+\label{introFig:eulerAngles}
+\end{figure}
+
+The equations of motion for Euler angles can be written down
+as\cite{allen87:csl}
+\begin{align}
+\dot{\phi} &= -\omega^s_x \frac{\sin\phi\cos\theta}{\sin\theta} +
+	\omega^s_y \frac{\cos\phi\cos\theta}{\sin\theta} +
+	\omega^s_z
+\label{introEq:MDeulerPhi} \\%
+%
+\dot{\theta} &= \omega^s_x \cos\phi + \omega^s_y \sin\phi
+\label{introEq:MDeulerTheta} \\%
+%
+\dot{\psi} &= \omega^s_x \frac{\sin\phi}{\sin\theta} - 
+	\omega^s_y \frac{\cos\phi}{\sin\theta}
+\label{introEq:MDeulerPsi}
+\end{align}
 Where $\omega^s_i$ is the angular velocity in the lab space frame
 along Cartesian coordinate $i$.  However, a difficulty arises when
 attempting to integrate Eq.~\ref{introEq:MDeulerPhi} and
 Eq.~\ref{introEq:MDeulerPsi}. The $\frac{1}{\sin \theta}$ present in
 both equations means there is a non-physical instability present when
-$\theta$ is 0 or $\pi$.
-
-To correct for this, the simulations integrate the rotation matrix,
-$\mathbf{A}$, directly, thus avoiding the instability.
-This method was proposed by Dullwebber
-\emph{et. al.}\cite{Dullwebber:1997}, and is presented in
+$\theta$ is 0 or $\pi$. To correct for this, the simulations integrate
+the rotation matrix, $\mathbf{A}$, directly, thus avoiding the
+instability.  This method was proposed by Dullweber
+\emph{et. al.}\cite{Dullweber1997}, and is presented in
 Sec.~\ref{introSec:MDsymplecticRot}.
 
-\subsubsection{\label{introSec:MDliouville}Liouville Propagator}
+\subsection{\label{introSec:MDliouville}Liouville Propagator}
 
 Before discussing the integration of the rotation matrix, it is
 necessary to understand the construction of a ``good'' integration
 scheme.  It has been previously
-discussed(Sec.~\ref{introSec:MDintegrate}) how it is desirable for an
+discussed(Sec.~\ref{introSec:mdIntegrate}) how it is desirable for an
 integrator to be symplectic, or time reversible.  The following is an
 outline of the Trotter factorization of the Liouville Propagator as a
-scheme for generating symplectic integrators. \cite{Tuckerman:1992}
+scheme for generating symplectic integrators.\cite{Tuckerman92}
 
 For a system with $f$ degrees of freedom the Liouville operator can be
 defined as,
 \begin{equation}
-eq here
+iL=\sum^f_{j=1} \biggl [\dot{q}_j\frac{\partial}{\partial q_j} +
+	F_j\frac{\partial}{\partial p_j} \biggr ]
 \label{introEq:LiouvilleOperator}
 \end{equation}
-Here, $r_j$ and $p_j$ are the position and conjugate momenta of a
-degree of freedom, and $f_j$ is the force on that degree of freedom.
+Here, $q_j$ and $p_j$ are the position and conjugate momenta of a
+degree of freedom, and $F_j$ is the force on that degree of freedom.
 $\Gamma$ is defined as the set of all positions and conjugate momenta,
-$\{r_j,p_j\}$, and the propagator, $U(t)$, is defined
+$\{q_j,p_j\}$, and the propagator, $U(t)$, is defined
 \begin {equation}
-eq here
+U(t) = e^{iLt}
 \label{introEq:Lpropagator}
 \end{equation}
 This allows the specification of $\Gamma$ at any time $t$ as 
 \begin{equation}
-eq here
+\Gamma(t) = U(t)\Gamma(0)
 \label{introEq:Lp2}
 \end{equation}
 It is important to note, $U(t)$ is a unitary operator meaning
@@ -680,42 +712,72 @@ Trotter theorem to yield
 
 Decomposing $L$ into two parts, $iL_1$ and $iL_2$, one can use the
 Trotter theorem to yield
-\begin{equation}
-eq here
-\label{introEq:Lp4}
-\end{equation}
-Where $\Delta t = \frac{t}{P}$.
+\begin{align}
+e^{iLt} &= e^{i(L_1 + L_2)t} \notag \\%
+%
+	&= \biggl [ e^{i(L_1 +L_2)\frac{t}{P}} \biggr]^P \notag \\%
+%
+	&= \biggl [ e^{iL_1\frac{\Delta t}{2}}\, e^{iL_2\Delta t}\,
+	e^{iL_1\frac{\Delta t}{2}} \biggr ]^P +
+	\mathcal{O}\biggl (\frac{t^3}{P^2} \biggr ) \label{introEq:Lp4}
+\end{align}
+Where $\Delta t = t/P$.
 With this, a discrete time operator $G(\Delta t)$ can be defined:
-\begin{equation}
-eq here
+\begin{align}
+G(\Delta t) &= e^{iL_1\frac{\Delta t}{2}}\, e^{iL_2\Delta t}\,
+	e^{iL_1\frac{\Delta t}{2}} \notag \\%
+%
+	&= U_1 \biggl ( \frac{\Delta t}{2} \biggr )\, U_2 ( \Delta t )\,
+	U_1 \biggl ( \frac{\Delta t}{2} \biggr )
 \label{introEq:Lp5}
-\end{equation}
-Because $U_1(t)$ and $U_2(t)$ are unitary, $G|\Delta t)$ is also
+\end{align}
+Because $U_1(t)$ and $U_2(t)$ are unitary, $G(\Delta t)$ is also
 unitary.  Meaning an integrator based on this factorization will be
 reversible in time.
 
 As an example, consider the following decomposition of $L$:
+\begin{align}
+iL_1 &= \dot{q}\frac{\partial}{\partial q}%
+\label{introEq:Lp6a} \\%
+%
+iL_2 &= F(q)\frac{\partial}{\partial p}%
+\label{introEq:Lp6b}
+\end{align}
+This leads to propagator $G( \Delta t )$ as,
 \begin{equation}
-eq here
-\label{introEq:Lp6}
+G(\Delta t) =  e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}} \,
+	e^{\Delta t\,\dot{q}\frac{\partial}{\partial q}} \,
+	e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}}
+\label{introEq:Lp7}
 \end{equation}
-Operating $G(\Delta t)$ on $\Gamma)0)$, and utilizing the operator property
+Operating $G(\Delta t)$ on $\Gamma(0)$, and utilizing the operator property
 \begin{equation}
-eq here
+e^{c\frac{\partial}{\partial x}}\, f(x) = f(x+c)
 \label{introEq:Lp8}
 \end{equation}
-Where $c$ is independent of $q$.  One obtains the following:  
-\begin{equation}
-eq here
-\label{introEq:Lp8}
-\end{equation}
+Where $c$ is independent of $x$.  One obtains the following:  
+\begin{align}
+\dot{q}\biggl (\frac{\Delta t}{2}\biggr ) &=
+	\dot{q}(0) + \frac{\Delta t}{2m}\, F[q(0)] \label{introEq:Lp9a}\\%
+%
+q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{\Delta t}{2}\biggr )%
+	\label{introEq:Lp9b}\\%
+%
+\dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
+	\frac{\Delta t}{2m}\, F[q(0)] \label{introEq:Lp9c}
+\end{align}
 Or written another way,
-\begin{equation}
-eq here
-\label{intorEq:Lp9}
-\end{equation}
+\begin{align}
+q(t+\Delta t) &= q(0) + \dot{q}(0)\Delta t + 
+	\frac{F[q(0)]}{m}\frac{\Delta t^2}{2} %
+\label{introEq:Lp10a} \\%
+%
+\dot{q}(\Delta t) &= \dot{q}(0) + \frac{\Delta t}{2m}
+	\biggl [F[q(0)] + F[q(\Delta t)] \biggr] %
+\label{introEq:Lp10b}
+\end{align}
 This is the velocity Verlet formulation presented in
-Sec.~\ref{introSec:MDintegrate}.  Because this integration scheme is
+Sec.~\ref{introSec:mdIntegrate}.  Because this integration scheme is
 comprised of unitary propagators, it is symplectic, and therefore area
 preserving in phase space.  From the preceding factorization, one can
 see that the integration of the equations of motion would follow:
@@ -729,70 +791,86 @@ see that the integration of the equations of motion wo
 \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
 \end{enumerate}
 
-\subsubsection{\label{introSec:MDsymplecticRot} Symplectic Propagation of the Rotation Matrix}
+\subsection{\label{introSec:MDsymplecticRot} Symplectic Propagation of the Rotation Matrix}
 
 Based on the factorization from the previous section,
-Dullweber\emph{et al.}\cite{Dullweber:1997}~ proposed a scheme for the
+Dullweber\emph{et al}.\cite{Dullweber1997}~ proposed a scheme for the
 symplectic propagation of the rotation matrix, $\mathbf{A}$, as an
 alternative method for the integration of orientational degrees of
 freedom. The method starts with a straightforward splitting of the
 Liouville operator:
-\begin{equation}
-eq here
-\label{introEq:SR1}
-\end{equation}
-Where $\boldsymbol{\tau}(\mathbf{A})$ are the torques of the system
-due to the configuration, and $\boldsymbol{/pi}$ are the conjugate
+\begin{align}
+iL_{\text{pos}} &= \dot{q}\frac{\partial}{\partial q} +
+	\mathbf{\dot{A}}\frac{\partial}{\partial \mathbf{A}} 
+\label{introEq:SR1a} \\%
+%
+iL_F &= F(q)\frac{\partial}{\partial p}
+\label{introEq:SR1b} \\%
+iL_{\tau} &= \tau(\mathbf{A})\frac{\partial}{\partial \pi}
+\label{introEq:SR1b} \\%
+\end{align}
+Where $\tau(\mathbf{A})$ is the torque of the system
+due to the configuration, and $\pi$ is the conjugate
 angular momenta of the system. The propagator, $G(\Delta t)$, becomes
 \begin{equation}
-eq here
+G(\Delta t) = e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}} \,
+	e^{\frac{\Delta t}{2} \tau(\mathbf{A})\frac{\partial}{\partial \pi}} \,
+	e^{\Delta t\,iL_{\text{pos}}} \,
+	e^{\frac{\Delta t}{2} \tau(\mathbf{A})\frac{\partial}{\partial \pi}} \,
+	e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}}
 \label{introEq:SR2}
 \end{equation}
 Propagation of the linear and angular momenta follows as in the Verlet
 scheme.  The propagation of positions also follows the Verlet scheme
 with the addition of a further symplectic splitting of the rotation
-matrix propagation, $\mathcal{G}_{\text{rot}}(\Delta t)$.
+matrix propagation, $\mathcal{U}_{\text{rot}}(\Delta t)$, within
+$U_{\text{pos}}(\Delta t)$.
 \begin{equation}
-eq here
+\mathcal{U}_{\text{rot}}(\Delta t) = 
+	\mathcal{U}_x \biggl(\frac{\Delta t}{2}\biggr)\,
+	\mathcal{U}_y \biggl(\frac{\Delta t}{2}\biggr)\,
+	\mathcal{U}_z (\Delta t)\,
+	\mathcal{U}_y \biggl(\frac{\Delta t}{2}\biggr)\,
+	\mathcal{U}_x \biggl(\frac{\Delta t}{2}\biggr)\,
 \label{introEq:SR3}
 \end{equation}
-Where $\mathcal{G}_j$ is a unitary rotation of $\mathbf{A}$ and
-$\boldsymbol{\pi}$ about each axis $j$.  As all propagations are now
+Where $\mathcal{U}_j$ is a unitary rotation of $\mathbf{A}$ and
+$\pi$ about each axis $j$.  As all propagations are now
 unitary and symplectic, the entire integration scheme is also
 symplectic and time reversible.
 
 \section{\label{introSec:layout}Dissertation Layout}
 
-This dissertation is divided as follows:Chapt.~\ref{chapt:RSA}
+This dissertation is divided as follows:Ch.~\ref{chapt:RSA}
 presents the random sequential adsorption simulations of related
 pthalocyanines on a gold (111) surface. Ch.~\ref{chapt:OOPSE}
 is about the writing of the molecular dynamics simulation package
-{\sc oopse}, Ch.~\ref{chapt:lipid} regards the simulations of
-phospholipid bilayers using a mesoscale model, and lastly,
+{\sc oopse}. Ch.~\ref{chapt:lipid} regards the simulations of
+phospholipid bilayers using a mesoscale model. And lastly,
 Ch.~\ref{chapt:conclusion} concludes this dissertation with a
 summary of all results. The chapters are arranged in chronological
 order, and reflect the progression of techniques I employed during my
 research.  
 
-The chapter concerning random sequential adsorption
-simulations is a study in applying the principles of theoretical
-research in order to obtain a simple model capable of explaining the
-results.  My advisor, Dr. Gezelter, and I were approached by a
-colleague, Dr. Lieberman, about possible explanations for partial
-coverage of a gold surface by a particular compound of hers. We
-suggested it might be due to the statistical packing fraction of disks
-on a plane, and set about to simulate this system.  As the events in
-our model were not dynamic in nature, a Monte Carlo method was
-employed.  Here, if a molecule landed on the surface without
-overlapping another, then its landing was accepted.  However, if there
-was overlap, the landing we rejected and a new random landing location
-was chosen.  This defined our acceptance rules and allowed us to
-construct a Markov chain whose limiting distribution was the surface
-coverage in which we were interested.
+The chapter concerning random sequential adsorption simulations is a
+study in applying Statistical Mechanics simulation techniques in order
+to obtain a simple model capable of explaining the results.  My
+advisor, Dr. Gezelter, and I were approached by a colleague,
+Dr. Lieberman, about possible explanations for the  partial coverage of a
+gold surface by a particular compound of hers. We suggested it might
+be due to the statistical packing fraction of disks on a plane, and
+set about to simulate this system.  As the events in our model were
+not dynamic in nature, a Monte Carlo method was employed.  Here, if a
+molecule landed on the surface without overlapping another, then its
+landing was accepted.  However, if there was overlap, the landing we
+rejected and a new random landing location was chosen.  This defined
+our acceptance rules and allowed us to construct a Markov chain whose
+limiting distribution was the surface coverage in which we were
+interested.
 
 The following chapter, about the simulation package {\sc oopse},
 describes in detail the large body of scientific code that had to be
-written in order to study phospholipid bilayer.  Although there are
+written in order to study phospholipid bilayers.  Although there are
 pre-existing molecular dynamic simulation packages available, none
 were capable of implementing the models we were developing.{\sc oopse}
 is a unique package capable of not only integrating the equations of
@@ -804,9 +882,9 @@ able to parameterize a mesoscale model for phospholipi
 
 Bringing us to Ch.~\ref{chapt:lipid}. Using {\sc oopse}, I have been
 able to parameterize a mesoscale model for phospholipid simulations.
-This model retains information about solvent ordering about the
+This model retains information about solvent ordering around the
 bilayer, as well as information regarding the interaction of the
-phospholipid head groups' dipole with each other and the surrounding
+phospholipid head groups' dipoles with each other and the surrounding
 solvent.  These simulations give us insight into the dynamic events
 that lead to the formation of phospholipid bilayers, as well as
 provide the foundation for future exploration of bilayer phase
