[ViewVC] Diff of: group/trunk/mattDisertation/Introduction.tex

Comparing trunk/mattDisertation/Introduction.tex (file contents):
Revision 977 by mmeineke, Thu Jan 22 21:13:55 2004 UTC vs.
Revision 1001 by mmeineke, Sat Jan 31 22:10:21 2004 UTC

+thermodynamic properties of the system are being probed, then chose
+which method best suits that objective.
-<
+\subsection{\label{introSec:statThermo}Statistical Thermodynamics}
->
+\subsection{\label{introSec:statThermo}Statistical Mechanics}
-<
+ergodic hypothesis
->
+The following section serves as a brief introduction to some of the
->
+Statistical Mechanics concepts present in this dissertation.  What
->
+follows is a brief derivation of Blotzman weighted statistics, and an
->
+explanation of how one can use the information to calculate an
->
+observable for a system.  This section then concludes with a brief
->
+discussion of the ergodic hypothesis and its relevance to this
->
+research.
-<
+enesemble averages
->
+\subsection{\label{introSec:boltzman}Boltzman weighted statistics}
-+
+Consider a system, $\gamma$, with some total energy,, $E_{\gamma}$.
-+
+Let $\Omega(E_{gamma})$ represent the number of degenerate ways
-+
+$\boldsymbol{\Gamma}$, the collection of positions and conjugate
-+
+momenta of system $\gamma$, can be configured to give
-+
+$E_{\gamma}$. Further, if $\gamma$ is in contact with a bath system
-+
+where energy is exchanged between the two systems, $\Omega(E)$, where
-+
+$E$ is the total energy of both systems, can be represented as
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM1}
-+
+\end{equation}
-+
+Or additively as
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM2}
-+
+\end{equation}
-+
-+
+The solution to Eq.~\ref{introEq:SM2} maximizes the number of
-+
+degenerative configurations in $E$. \cite{fix}
-+
+This gives
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM3}
-+
+\end{equation}
-+
+Where $E_{\text{bath}}$ is $E-E_{\gamma}$, and
-+
+$\frac{partialE_{\text{bath}}}{\partial E_{\gamma}}$ is
-+
+$-1$. Eq.~\ref{introEq:SM3} becomes
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM4}
-+
+\end{equation}
-+
-+
+At this point, one can draw a relationship between the maximization of
-+
+degeneracy in Eq.~\ref{introEq:SM3} and the second law of
-+
+thermodynamics.  Namely, that for a closed system, entropy wil
-+
+increase for an irreversible process.\cite{fix} Here the
-+
+process is the partitioning of energy among the two systems.  This
-+
+allows the following definition of entropy:
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM5}
-+
+\end{equation}
-+
+Where $k_B$ is the Boltzman constant.  Having defined entropy, one can
-+
+also define the temperature of the system using the relation
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM6}
-+
+\end{equation}
-+
+The temperature in the system $\gamma$ is then
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM7}
-+
+\end{equation}
-+
+Applying this to Eq.~\ref{introEq:SM4} gives the following
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM8}
-+
+\end{equation}
-+
+Showing that the partitioning of energy between the two systems is
-+
+actually a process of thermal equilibration. \cite{fix}
-+
-+
+An application of these results is to formulate the form of an
-+
+expectation value of an observable, $A$, in the canonical ensemble. In
-+
+the canonical ensemble, the number of particles, $N$, the volume, $V$,
-+
+and the temperature, $T$, are all held constant while the energy, $E$,
-+
+is allowed to fluctuate. Returning to the previous example, the bath
-+
+system is now an infinitly large thermal bath, whose exchange of
-+
+energy with the system $\gamma$ holds teh temperature constant.  The
-+
+partitioning of energy in the bath system is then related to the total
-+
+energy of both systems and the fluctuations in $E_{\gamma}}$:
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM9}
-+
+\end{equation}
-+
+As for the expectation value, it can be defined
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM10}
-+
+\end{eequation}
-+
+Where $\int_{\boldsymbol{\Gamma}} d\Boldsymbol{\Gamma}$ denotes an
-+
+integration over all accessable phase space, $P_{\gamma}$ is the
-+
+probability of being in a given phase state and
-+
+$A(\boldsymbol{\Gamma})$ is some observable that is a function of the
-+
+phase state.
-+
-+
+Because entropy seeks to maximize the number of degenerate states at a
-+
+given energy, the probability of being in a particular state in
-+
+$\gamma$ will be directly proportional to the number of allowable
-+
+states the coupled system is able to assume. Namely,
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM11}
-+
+\end{equation}
-+
+With $E_{\gamma} \lE$, $\ln \Omega$ can be expanded in a Taylor series:
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM12}
-+
+\end{equation}
-+
+Higher order terms are omitted as $E$ is an infinite thermal
-+
+bath. Further, using Eq.~\ref{introEq:SM7}, Eq.~\ref{introEq:SM11} can
-+
+be rewritten:
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM13}
-+
+\end{equation}
-+
+Where $\ln \Omega(E)$ has been factored out of the porpotionality as a
-+
+constant.  Normalizing the probability ($\int_{\boldsymbol{\Gamma}}
-+
+d\boldsymbol{\Gamma} P_{\gamma} =1$ gives
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM14}
-+
+\end{equation}
-+
+This result is the standard Boltzman statistical distribution.
-+
+Applying it to Eq.~\ref{introEq:SM10} one can obtain the following relation for ensemble averages:
-+
+\begin{equation}
-+
+eq here
-+
+\label{introEq:SM15}
-+
+\end{equation}
-+
-+
+\subsection{\label{introSec:ergodic}The Ergodic Hypothesis}
-+
-+
+One last important consideration is that of ergodicity. Ergodicity is
-+
+the assumption that given an infinite amount of time, a system will
-+
+visit every available point in phase space.\cite{fix} For most
-+
+systems, this is a valid assumption, except in cases where the system
-+
+may be trapped in a local feature (\emph{i.~e.~glasses}). When valid,
-+
+ergodicity allows the unification of a time averaged observation and
-+
+an ensemble averged one. If an observation is averaged over a
-+
+sufficiently long time, the system is assumed to visit all
-+
+appropriately available points in phase space, giving a properly
-+
+weighted statistical average. This allows the researcher freedom of
-+
+choice when deciding how best to measure a given observable.  When an
-+
+ensemble averaged approach seems most logical, the Monte Carlo
-+
+techniques described in Sec.~\ref{introSec:MC} can be utilized.
-+
+Conversely, if a problem lends itself to a time averaging approach,
-+
+the Molecular Dynamics techniques in Sec.~\ref{introSec:MD} can be
-+
+employed.
-+
+\subsection{\label{introSec:monteCarlo}Monte Carlo Simulations}
+The Monte Carlo method was developed by Metropolis and Ulam for their
+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
->
+researcher is interested.  If the observables 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.
+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
->
+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}
->
+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
+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$
->
+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}).
-+
+With the use of an $fix$, however, comes a discontinuity in the
-+
+potential energy curve (Fig.~\ref{fix}). To fix this discontinuity,
-+
+one calculates the potential energy at the $r_{\text{cut}}$, and add
-+
+that value to the potential.  This causes the function to go smoothly
-+
+to zero at the cutoff radius.  This ensures conservation of energy
-+
+when integrating the Newtonian equations of motion.
-<
+\section{\label{introSec:chapterLayout}Chapter Layout}
-<
-<
+\subsection{\label{introSec:RSA}Random Sequential Adsorption}
->
+The second main simplification used in this research is the Verlet
->
+neighbor list. \cite{allen87:csl} In the Verlet method, one generates
->
+a list of all neighbor atoms, $j$, surrounding atom $i$ within some
->
+cutoff $r_{\text{list}}$, where $r_{\text{list}}>r_{\text{cut}}$.
->
+This list is created the first time forces are evaluated, then on
->
+subsequent force evaluations, pair calculations are only calculated
->
+from the neighbor lists.  The lists are updated if any given particle
->
+in the system moves farther than $r_{\text{list}}-r_{\text{cut}}$,
->
+giving rise to the possibility that a particle has left or joined a
->
+neighbor list.
-<
+\subsection{\label{introSec:OOPSE}The OOPSE Simulation Package}
->
+\subsection{\label{introSec:MDintegrate} Integration of the equations of motion}
->
->
+A starting point for the discussion of molecular dynamics integrators
->
+is the Verlet algorithm. \cite{Frenkel1996} It begins with a Taylor
->
+expansion of position in time:
->
+\begin{equation}
->
+eq here
->
+\label{introEq:verletForward}
->
+\end{equation}
->
+As well as,
->
+\begin{equation}
->
+eq here
->
+\label{introEq:verletBack}
->
+\end{equation}
->
+Adding together Eq.~\ref{introEq:verletForward} and
->
+Eq.~\ref{introEq:verletBack} results in,
->
+\begin{equation}
->
+eq here
->
+\label{introEq:verletSum}
->
+\end{equation}
->
+Or equivalently,
->
+\begin{equation}
->
+eq here
->
+\label{introEq:verletFinal}
->
+\end{equation}
->
+Which contains an error in the estimate of the new positions on the
->
+order of $\Delta t^4$.
->
->
+In practice, however, the simulations in this research were integrated
->
+with a velocity reformulation of teh Verlet method.\cite{allen87:csl}
->
+\begin{equation}
->
+eq here
->
+\label{introEq:MDvelVerletPos}
->
+\end{equation}
->
+\begin{equation}
->
+eq here
->
+\label{introEq:MDvelVerletVel}
->
+\end{equation}
->
+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
->
+very little long term drift in energy conservation.  Energy
->
+conservation in a molecular dynamics simulation is of extreme
->
+importance, as it is a measure of how closely one is following the
->
+``true'' trajectory wtih the finite integration scheme.  An exact
->
+solution to the integration will conserve area in phase space, as well
->
+as be reversible in time, that is, the trajectory integrated forward
->
+or backwards will exactly match itself.  Having a finite algorithm
->
+that both conserves area in phase space and is time reversible,
->
+therefore increases, but does not guarantee the ``correctness'' or the
->
+integrated trajectory.
-<
+\subsection{\label{introSec:bilayers}A Mesoscale Model for Phospholipid Bilayers}
->
+It can be shown,\cite{Frenkel1996} that although the Verlet algorithm
->
+does not rigorously preserve the actual Hamiltonian, it does preserve
->
+a pseudo-Hamiltonian which shadows the real one in phase space.  This
->
+pseudo-Hamiltonian is proveably area-conserving as well as time
->
+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
->
+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 Chapt.~\ref{chaptLipids} 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}$
->
+defined as follows:
->
+\begin{equation}
->
+eq here
->
+\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}
->
+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
->
+Sec.~\ref{introSec:MDsymplecticRot}.
->
->
+\subsubsection{\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
->
+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}
->
->
+For a system with $f$ degrees of freedom the Liouville operator can be
->
+defined as,
->
+\begin{equation}
->
+eq here
->
+\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.
->
+$\Gamma$ is defined as the set of all positions nad conjugate momenta,
->
+$\{r_j,p_j\}$, and the propagator, $U(t)$, is defined
->
+\begin {equation}
->
+eq here
->
+\label{introEq:Lpropagator}
->
+\end{equation}
->
+This allows the specification of $\Gamma$ at any time $t$ as
->
+\begin{equation}
->
+eq here
->
+\label{introEq:Lp2}
->
+\end{equation}
->
+It is important to note, $U(t)$ is a unitary operator meaning
->
+\begin{equation}
->
+U(-t)=U^{-1}(t)
->
+\label{introEq:Lp3}
->
+\end{equation}
->
->
+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}$.
->
+With this, a discrete time operator $G(\Delta t)$ can be defined:
->
+\begin{equation}
->
+eq here
->
+\label{introEq:Lp5}
->
+\end{equation}
->
+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{equation}
->
+eq here
->
+\label{introEq:Lp6}
->
+\end{equation}
->
+Operating $G(\Delta t)$ on $\Gamma)0)$, and utilizing the operator property
->
+\begin{equation}
->
+eq here
->
+\label{introEq:Lp8}
->
+\end{equation}
->
+Where $c$ is independent of $q$.  One obtains the following:
->
+\begin{equation}
->
+eq here
->
+\label{introEq:Lp8}
->
+\end{equation}
->
+Or written another way,
->
+\begin{equation}
->
+eq here
->
+\label{intorEq:Lp9}
->
+\end{equation}
->
+This is the velocity Verlet formulation presented in
->
+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 preceeding fatorization, one can
->
+see that the integration of the equations of motion would follow:
->
+\begin{enumerate}
->
+\item calculate the velocities at the half step, $\frac{\Delta t}{2}$, from the forces calculated at the initial position.
->
->
+\item Use the half step velocities to move positions one whole step, $\Delta t$.
->
->
+\item Evaluate the forces at the new positions, $\mathbf{r}(\Delta t)$, and use the new forces to complete the velocity move.
->
->
+\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}
->
->
+Based on the factorization from the previous section,
->
+Dullweber\emph{et al.}\cite{Dullweber:1997}~ 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 tourques of the system
->
+due to the configuration, and $\boldsymbol{/pi}$ are the conjugate
->
+angular momenta of the system. The propagator, $G(\Delta t)$, becomes
->
+\begin{equation}
->
+eq here
->
+\label{introEq:SR2}
->
+\end{equation}
->
+Propagation fo 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)$.
->
+\begin{equation}
->
+eq here
->
+\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
->
+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}
->
+presents the random sequential adsorption simulations of related
->
+pthalocyanines on a gold (111) surface. Chapt.~\ref{chapt:OOPSE}
->
+is about the writing of the molecular dynamics simulation package
->
+{\sc oopse}, Chapt.~\ref{chapt:lipid} regards the simulations of
->
+phospholipid bilayers using a mesoscale model, and lastly,
->
+Chapt.~\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 capaable of explaining the
->
+results.  My advisor, Dr. Gezelter, and I were approached by a
->
+colleague, Dr. Lieberman, about possible explanations for partial
->
+coverge 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
->
+emplyed.  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
->
+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
->
+motion in cartesian space, but is also able to integrate the
->
+rotational motion of rigid bodies and dipoles.  Add to this the
->
+ability to perform calculations across parallel processors and a
->
+flexible script syntax for creating systems, and {\sc oopse} becomes a
->
+very powerful scientific instrument for the exploration of our model.
->
->
+Bringing us to Chapt.~\ref{chapt:lipid}. Using {\sc oopse}, I have been
->
+able to parametrize a mesoscale model for phospholipid simulations.
->
+This model retains information about solvent ordering about the
->
+bilayer, as well as information regarding the interaction of the
->
+phospholipid head groups' dipole 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
->
+behavior with this model.
->
->
+Which leads into the last chapter, where I discuss future directions
->
+for both{\sc oopse} and this mesoscale model.  Additionally, I will
->
+give a summary of results for this dissertation.
->
->

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing trunk/mattDisertation/Introduction.tex (file contents): Revision 977 by mmeineke, Thu Jan 22 21:13:55 2004 UTC vs. Revision 1001 by mmeineke, Sat Jan 31 22:10:21 2004 UTC

Diff Legend

Comparing trunk/mattDisertation/Introduction.tex (file contents):
Revision 977 by mmeineke, Thu Jan 22 21:13:55 2004 UTC vs.
Revision 1001 by mmeineke, Sat Jan 31 22:10:21 2004 UTC