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

Comparing trunk/mattDisertation/Introduction.tex (file contents):
Revision 953 by mmeineke, Sun Jan 18 03:07:59 2004 UTC vs.
Revision 1003 by mmeineke, Mon Feb 2 21:56:16 2004 UTC

+\chapter{\label{chapt:intro}Introduction and Theoretical Background}
-–
-–
+\section{\label{introSec:theory}Theoretical Background}
-–
+The techniques used in the course of this research fall under the two
+main classes of molecular simulation: Molecular Dynamics and Monte
+Carlo. Molecular Dynamic simulations integrate the equations of motion
+thermodynamic properties of the system are being probed, then chose
+which method best suits that objective.
-<
+\subsection{\label{introSec:statThermo}Statistical Thermodynamics}
->
+\section{\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}
-<
+\subsection{\label{introSec:monteCarlo}Monte Carlo Simulations}
->
+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}
->
+\Omega(E) = \Omega(E_{\gamma}) \times \Omega(E - E_{\gamma})
->
+\label{introEq:SM1}
->
+\end{equation}
->
+Or additively as
->
+\begin{equation}
->
+\ln \Omega(E) = \ln \Omega(E_{\gamma}) + \ln \Omega(E - E_{\gamma})
->
+\label{introEq:SM2}
->
+\end{equation}
-+
+The solution to Eq.~\ref{introEq:SM2} maximizes the number of
-+
+degenerative configurations in $E$. \cite{Frenkel1996}
-+
+This gives
-+
+\begin{equation}
-+
+\frac{\partial \ln \Omega(E)}{\partial E_{\gamma}} = 0 =
-+
+        \frac{\partial \ln \Omega(E_{\gamma})}{\partial E_{\gamma}}
-+
+         +
-+
+        \frac{\partial \ln \Omega(E_{\text{bath}})}{\partial E_{\text{bath}}}
-+
+        \frac{\partial E_{\text{bath}}}{\partial E_{\gamma}}
-+
+\label{introEq:SM3}
-+
+\end{equation}
-+
+Where $E_{\text{bath}}$ is $E-E_{\gamma}$, and
-+
+$\frac{\partial E_{\text{bath}}}{\partial E_{\gamma}}$ is
-+
+$-1$. Eq.~\ref{introEq:SM3} becomes
-+
+\begin{equation}
-+
+\frac{\partial \ln \Omega(E_{\gamma})}{\partial E_{\gamma}} =
-+
+\frac{\partial \ln \Omega(E_{\text{bath}})}{\partial E_{\text{bath}}}
-+
+\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 will
-+
+increase for an irreversible process.\cite{chandler:1987} Here the
-+
+process is the partitioning of energy among the two systems.  This
-+
+allows the following definition of entropy:
-+
+\begin{equation}
-+
+S = k_B \ln \Omega(E)
-+
+\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}
-+
+\frac{1}{T} = \biggl ( \frac{\partial S}{\partial E} \biggr )_{N,V}
-+
+\label{introEq:SM6}
-+
+\end{equation}
-+
+The temperature in the system $\gamma$ is then
-+
+\begin{equation}
-+
+\beta( E_{\gamma} ) = \frac{1}{k_B T} =
-+
+        \frac{\partial \ln \Omega(E_{\gamma})}{\partial E_{\gamma}}
-+
+\label{introEq:SM7}
-+
+\end{equation}
-+
+Applying this to Eq.~\ref{introEq:SM4} gives the following
-+
+\begin{equation}
-+
+\beta( E_{\gamma} ) = \beta( E_{\text{bath}} )
-+
+\label{introEq:SM8}
-+
+\end{equation}
-+
+Showing that the partitioning of energy between the two systems is
-+
+actually a process of thermal equilibration.\cite{Frenkel1996}
-+
-+
+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 the 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}
-+
+\Omega( E_{\gamma} ) = \Omega( E - E_{\gamma} )
-+
+\label{introEq:SM9}
-+
+\end{equation}
-+
+As for the expectation value, it can be defined
-+
+\begin{equation}
-+
+\langle A \rangle =
-+
+        \int\limits_{\boldsymbol{\Gamma}} d\boldsymbol{\Gamma}
-+
+        P_{\gamma} A(\boldsymbol{\Gamma})
-+
+\label{introEq:SM10}
-+
+\end{equation}
-+
+Where $\int\limits_{\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}
-+
+P_{\gamma} \propto \Omega( E_{\text{bath}} ) =
-+
+        e^{\ln \Omega( E - E_{\gamma})}
-+
+\label{introEq:SM11}
-+
+\end{equation}
-+
+With $E_{\gamma} \ll E$, $\ln \Omega$ can be expanded in a Taylor series:
-+
+\begin{equation}
-+
+\ln \Omega ( E - E_{\gamma}) =
-+
+        \ln \Omega (E) -
-+
+        E_{\gamma}  \frac{\partial \ln \Omega }{\partial E}
-+
+        + \ldots
-+
+\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}
-+
+P_{\gamma} \propto e^{-\beta E_{\gamma}}
-+
+\label{introEq:SM13}
-+
+\end{equation}
-+
+Where $\ln \Omega(E)$ has been factored out of the porpotionality as a
-+
+constant.  Normalizing the probability ($\int\limits_{\boldsymbol{\Gamma}}
-+
+d\boldsymbol{\Gamma} P_{\gamma} = 1$) gives
-+
+\begin{equation}
-+
+P_{\gamma} = \frac{e^{-\beta E_{\gamma}}}
-+
+{\int\limits_{\boldsymbol{\Gamma}} d\boldsymbol{\Gamma} e^{-\beta E_{\gamma}}}
-+
+\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}
-+
+\langle A \rangle =
-+
+\frac{\int\limits_{\boldsymbol{\Gamma}} d\boldsymbol{\Gamma}
-+
+        A( \boldsymbol{\Gamma} ) e^{-\beta E_{\gamma}}}
-+
+{\int\limits_{\boldsymbol{\Gamma}} d\boldsymbol{\Gamma} e^{-\beta E_{\gamma}}}
-+
+\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{Frenkel1996} For most
-+
+systems, this is a valid assumption, except in cases where the system
-+
+may be trapped in a local feature (\emph{e.g.}~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:monteCarlo} 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.
-+
-+
+\section{\label{introSec:monteCarlo}Monte Carlo Simulations}
-+
+The Monte Carlo method was developed by Metropolis and Ulam for their
+work in fissionable material.\cite{metropolis:1949} The method is so
-<
+named, because it heavily uses random numbers in the solution of the
-<
+problem.
->
+named, because it heavily uses random numbers in its
->
+solution.\cite{allen87:csl} The Monte Carlo method allows for the
->
+solution of integrals through the stochastic sampling of the values
->
+within the integral. In the simplest case, the evaluation of an
->
+integral would follow a brute force method of
->
+sampling.\cite{Frenkel1996} Consider the following single dimensional
->
+integral:
->
+\begin{equation}
->
+I = f(x)dx
->
+\label{eq:MCex1}
->
+\end{equation}
->
+The equation can be recast as:
->
+\begin{equation}
->
+I = (b-a)\langle f(x) \rangle
->
+\label{eq:MCex2}
->
+\end{equation}
->
+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
->
+approach $I$ in the limit where the number of trials is infintely
->
+large.
-+
+However, in Statistical Mechanics, one is typically interested in
-+
+integrals of the form:
-+
+\begin{equation}
-+
+\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 $\mathbf{r}^N$ stands for the coordinates of all $N$ particles
-+
+and $A$ is some observable that is only dependent on position. This is
-+
+the ensemble average of $A$ as presented in
-+
+Sec.~\ref{introSec:statThermo}, except here $A$ is independent of
-+
+momentum. Therefore 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 importance sampling comes into
-+
+play.\cite{allen87:csl}
-<
+\subsection{\label{introSec:md}Molecular Dynamics Simulations}
->
+Importance Sampling is a method where one selects a distribution from
->
+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)$.
-<
+time averages
->
+\subsection{\label{introSec:markovChains}Markov Chains}
-<
+time integrating schemes
->
+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$.
-<
+time reversible
->
+\newcommand{\accMe}{\operatorname{acc}}
-<
+symplectic methods
->
+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}
-<
+Extended ensembles (NVT NPT)
->
+\subsection{\label{introSec:metropolisMethod}The Metropolis Method}
-<
+constrained dynamics
->
+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}
->
+\accMe( m \rightarrow n ) =
->
+        \begin{cases}
->
+& \text{if $\Delta \mathcal{U} \leq 0$,} \\
->
+        e^{-\beta \Delta \mathcal{U}}& \text{if $\Delta \mathcal{U} > 0$.}
->
+        \end{cases}
->
+\label{introEq:accRules}
->
+\end{equation}
-<
+\section{\label{introSec:chapterLayout}Chapter Layout}
->
+Using the acceptance criteria from Eq.~\ref{introEq:accRules} the
->
+Metropolis method proceeds as follows
->
+\begin{enumerate}
->
+\item Generate an initial configuration $\mathbf{r}^N$ which has some finite probability in $\rho_{kT}$.
->
+\item Modify $\mathbf{r}^N$, to generate configuratioon $\mathbf{r^{\prime}}^N$.
->
+\item If the new configuration lowers the energy of the system, accept the move with unity ($\mathbf{r}^N$ becomes $\mathbf{r^{\prime}}^N$).  Otherwise accept with probability $e^{-\beta \Delta \mathcal{U}}$.
->
+\item Accumulate the average for the configurational observable of intereest.
->
+\item Repeat from step 2 until the average converges.
->
+\end{enumerate}
->
+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{introEq:Importance4} it becomes clear that the
->
+accumulated averages are the ensemble averages, as this method ensures
->
+that the limiting distribution is the Boltzman distribution.
-<
+\subsection{\label{introSec:RSA}Random Sequential Adsorption}
->
+\section{\label{introSec:MD}Molecular Dynamics Simulations}
-<
+\subsection{\label{introSec:OOPSE}The OOPSE Simulation Package}
->
+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,
->
+Sec.~\ref{introSec:ergodic}, the average of these observables over the
->
+time period of the simulation are taken to be the ensemble averages
->
+for the system.
-<
+\subsection{\label{introSec:bilayers}A Mesoscale Model for Phospholipid Bilayers}
->
+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 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.
->
->
+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.
->
->
+\subsection{\label{introSec:mdAlgorithm}The 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{introSec:mdInit} deals with the initialization
->
+of a simulation.  Sec.~\ref{introSec:mdForce} discusses issues
->
+involved with the calculation of the forces.
->
+Sec.~\ref{introSec:mdIntegrate} concludes the algorithm discussion
->
+with the integration of the equations of motion.\cite{Frenkel1996}
->
->
+\subsection{\label{introSec:mdInit}Initialization}
->
->
+When selecting the initial configuration for the simulation it is
->
+important to consider what dynamics one is hoping to observe.
->
+Ch.~\ref{chapt:lipid} 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 the 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.
->
->
+\subsection{\label{introSec:mdForce}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 $N(N-1)/ 2$
->
+pairs to be evaluated, where $N$ is the number of particles in the
->
+system.  This leads to the calculations scaling as $N^2$, 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{allen87:csl} For a cubic system of 1000 particles
->
+arranged in a $10 \times 10 \times 10$ 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{born:1912}
->
->
+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{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
->
+Sec.\ref{oopseSec:pbc}.
->
->
+\begin{figure}
->
+\centering
->
+\includegraphics[width=\linewidth]{pbcFig.eps}
->
+\caption[An illustration of periodic boundry conditions]{A 2-D illustration of periodic boundry conditions. As one particle leaves the right of the simulation box, an image of it enters the left.}
->
+\label{introFig:pbc}
->
+\end{figure}
->
->
+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, $r_{\text{cut}}$, are not included in the
->
+calculation.\cite{Frenkel1996} In a simultation with periodic images,
->
+$r_{\text{cut}}$ has a maximum value of $\text{box}/2$.
->
+Fig.~\ref{introFig:rMax} illustrates how when using an
->
+$r_{\text{cut}}$ larger than this value, or in the extreme limit of no
->
+$r_{\text{cut}}$ 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 $\text{box} / 2$.
->
->
+\begin{figure}
->
+\centering
->
+\includegraphics[width=\linewidth]{rCutMaxFig.eps}
->
+\caption
->
+\label{introFig:rMax}
->
+\end{figure}
->
->
+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.
->
->
+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: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.
->
->
+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 953 by mmeineke, Sun Jan 18 03:07:59 2004 UTC vs. Revision 1003 by mmeineke, Mon Feb 2 21:56:16 2004 UTC

Diff Legend

Comparing trunk/mattDisertation/Introduction.tex (file contents):
Revision 953 by mmeineke, Sun Jan 18 03:07:59 2004 UTC vs.
Revision 1003 by mmeineke, Mon Feb 2 21:56:16 2004 UTC