ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mattDisertation/Introduction.tex
(Generate patch)

Comparing trunk/mattDisertation/Introduction.tex (file contents):
Revision 914 by mmeineke, Fri Jan 9 02:41:45 2004 UTC vs.
Revision 977 by mmeineke, Thu Jan 22 21:13:55 2004 UTC

# Line 6 | Line 6
6  
7   \section{\label{introSec:theory}Theoretical Background}
8  
9 < The techniques used in the course of this research fall under the two main classes of
10 < molecular simulation: Molecular Dynamics and Monte Carlo. Molecular Dynamic simulations
11 < integrate the equations of motion for a given system of particles, allowing the researher
12 < to gain insight into the time dependent evolution of a system. Diffusion phenomena are
13 < readily studied with this simulation technique, making Molecular Dynamics the main simulation
14 < technique used in this research. Other aspects of the research fall under the Monte Carlo
15 < class of simulations. In Monte Carlo, the configuration space available to the collection
16 < of particles is sampled stochastichally, or randomly. Each configuration is chosen with
17 < a given probability based on the Maxwell Boltzman distribution. These types of simulations
18 < are best used to probe properties of a system that are only dependent only on the state of
19 < the system. Structural information about a system is most readily obtained through
20 < these types of methods.
9 > The techniques used in the course of this research fall under the two
10 > main classes of molecular simulation: Molecular Dynamics and Monte
11 > Carlo. Molecular Dynamic simulations integrate the equations of motion
12 > for a given system of particles, allowing the researher to gain
13 > insight into the time dependent evolution of a system. Diffusion
14 > phenomena are readily studied with this simulation technique, making
15 > Molecular Dynamics the main simulation technique used in this
16 > research. Other aspects of the research fall under the Monte Carlo
17 > class of simulations. In Monte Carlo, the configuration space
18 > available to the collection of particles is sampled stochastichally,
19 > or randomly. Each configuration is chosen with a given probability
20 > based on the Maxwell Boltzman distribution. These types of simulations
21 > are best used to probe properties of a system that are only dependent
22 > only on the state of the system. Structural information about a system
23 > is most readily obtained through these types of methods.
24  
25 < Although the two techniques employed seem dissimilar, they are both linked by the overarching
26 < principles of Statistical Thermodynamics. Statistical Thermodynamics governs the behavior of
27 < both classes of simulations and dictates what each method can and cannot do. When investigating
28 < a system, one most first analyze what thermodynamic properties of the system are being probed,
29 < then chose which method best suits that objective.
25 > Although the two techniques employed seem dissimilar, they are both
26 > linked by the overarching principles of Statistical
27 > Thermodynamics. Statistical Thermodynamics governs the behavior of
28 > both classes of simulations and dictates what each method can and
29 > cannot do. When investigating a system, one most first analyze what
30 > thermodynamic properties of the system are being probed, then chose
31 > which method best suits that objective.
32  
33   \subsection{\label{introSec:statThermo}Statistical Thermodynamics}
34  
# Line 33 | Line 38 | enesemble averages
38  
39   \subsection{\label{introSec:monteCarlo}Monte Carlo Simulations}
40  
41 < Stochastic sampling
41 > The Monte Carlo method was developed by Metropolis and Ulam for their
42 > work in fissionable material.\cite{metropolis:1949} The method is so
43 > named, because it heavily uses random numbers in its
44 > solution.\cite{allen87:csl} The Monte Carlo method allows for the
45 > solution of integrals through the stochastic sampling of the values
46 > within the integral. In the simplest case, the evaluation of an
47 > integral would follow a brute force method of
48 > sampling.\cite{Frenkel1996} Consider the following single dimensional
49 > integral:
50 > \begin{equation}
51 > I = f(x)dx
52 > \label{eq:MCex1}
53 > \end{equation}
54 > The equation can be recast as:
55 > \begin{equation}
56 > I = (b-a)\langle f(x) \rangle
57 > \label{eq:MCex2}
58 > \end{equation}
59 > Where $\langle f(x) \rangle$ is the unweighted average over the interval
60 > $[a,b]$. The calculation of the integral could then be solved by
61 > randomly choosing points along the interval $[a,b]$ and calculating
62 > the value of $f(x)$ at each point. The accumulated average would then
63 > approach $I$ in the limit where the number of trials is infintely
64 > large.
65  
66 < detatiled balance
66 > However, in Statistical Mechanics, one is typically interested in
67 > integrals of the form:
68 > \begin{equation}
69 > \langle A \rangle = \frac{\int d^N \mathbf{r}~A(\mathbf{r}^N)%
70 >        e^{-\beta V(\mathbf{r}^N)}}%
71 >        {\int d^N \mathbf{r}~e^{-\beta V(\mathbf{r}^N)}}
72 > \label{eq:mcEnsAvg}
73 > \end{equation}
74 > Where $\mathbf{r}^N$ stands for the coordinates of all $N$ particles
75 > and $A$ is some observable that is only dependent on
76 > position. $\langle A \rangle$ is the ensemble average of $A$ as
77 > presented in Sec.~\ref{introSec:statThermo}. Because $A$ is
78 > independent of momentum, the momenta contribution of the integral can
79 > be factored out, leaving the configurational integral. Application of
80 > the brute force method to this system would yield highly inefficient
81 > results. Due to the Boltzman weighting of this integral, most random
82 > configurations will have a near zero contribution to the ensemble
83 > average. This is where a importance sampling comes into
84 > play.\cite{allen87:csl}
85  
86 < metropilis monte carlo
86 > Importance Sampling is a method where one selects a distribution from
87 > which the random configurations are chosen in order to more
88 > efficiently calculate the integral.\cite{Frenkel1996} Consider again
89 > Eq.~\ref{eq:MCex1} rewritten to be:
90 > \begin{equation}
91 > I = \int^b_a \frac{f(x)}{\rho(x)} \rho(x) dx
92 > \label{introEq:Importance1}
93 > \end{equation}
94 > Where $\rho(x)$ is an arbitrary probability distribution in $x$.  If
95 > one conducts $\tau$ trials selecting a random number, $\zeta_\tau$,
96 > from the distribution $\rho(x)$ on the interval $[a,b]$, then
97 > Eq.~\ref{introEq:Importance1} becomes
98 > \begin{equation}
99 > I= \biggl \langle \frac{f(x)}{\rho(x)} \biggr \rangle_{\text{trials}}
100 > \label{introEq:Importance2}
101 > \end{equation}
102 > Looking at Eq.~\ref{eq:mcEnsAvg}, and realizing
103 > \begin {equation}
104 > \rho_{kT}(\mathbf{r}^N) =
105 >        \frac{e^{-\beta V(\mathbf{r}^N)}}
106 >        {\int d^N \mathbf{r}~e^{-\beta V(\mathbf{r}^N)}}
107 > \label{introEq:MCboltzman}
108 > \end{equation}
109 > Where $\rho_{kT}$ is the boltzman distribution.  The ensemble average
110 > can be rewritten as
111 > \begin{equation}
112 > \langle A \rangle = \int d^N \mathbf{r}~A(\mathbf{r}^N)
113 >        \rho_{kT}(\mathbf{r}^N)
114 > \label{introEq:Importance3}
115 > \end{equation}
116 > Applying Eq.~\ref{introEq:Importance1} one obtains
117 > \begin{equation}
118 > \langle A \rangle = \biggl \langle
119 >        \frac{ A \rho_{kT}(\mathbf{r}^N) }
120 >        {\rho(\mathbf{r}^N)} \biggr \rangle_{\text{trials}}
121 > \label{introEq:Importance4}
122 > \end{equation}
123 > By selecting $\rho(\mathbf{r}^N)$ to be $\rho_{kT}(\mathbf{r}^N)$
124 > Eq.~\ref{introEq:Importance4} becomes
125 > \begin{equation}
126 > \langle A \rangle = \langle A(\mathbf{r}^N) \rangle_{\text{trials}}
127 > \label{introEq:Importance5}
128 > \end{equation}
129 > The difficulty is selecting points $\mathbf{r}^N$ such that they are
130 > sampled from the distribution $\rho_{kT}(\mathbf{r}^N)$.  A solution
131 > was proposed by Metropolis et al.\cite{metropolis:1953} which involved
132 > the use of a Markov chain whose limiting distribution was
133 > $\rho_{kT}(\mathbf{r}^N)$.
134  
135 < \subsection{\label{introSec:md}Molecular Dynamics Simulations}
135 > \subsubsection{\label{introSec:markovChains}Markov Chains}
136  
137 < time averages
137 > A Markov chain is a chain of states satisfying the following
138 > conditions:\cite{leach01:mm}
139 > \begin{enumerate}
140 > \item The outcome of each trial depends only on the outcome of the previous trial.
141 > \item Each trial belongs to a finite set of outcomes called the state space.
142 > \end{enumerate}
143 > If given two configuartions, $\mathbf{r}^N_m$ and $\mathbf{r}^N_n$,
144 > $\rho_m$ and $\rho_n$ are the probablilities of being in state
145 > $\mathbf{r}^N_m$ and $\mathbf{r}^N_n$ respectively.  Further, the two
146 > states are linked by a transition probability, $\pi_{mn}$, which is the
147 > probability of going from state $m$ to state $n$.
148  
149 < time integrating schemes
149 > \newcommand{\accMe}{\operatorname{acc}}
150  
151 < time reversible
151 > The transition probability is given by the following:
152 > \begin{equation}
153 > \pi_{mn} = \alpha_{mn} \times \accMe(m \rightarrow n)
154 > \label{introEq:MCpi}
155 > \end{equation}
156 > Where $\alpha_{mn}$ is the probability of attempting the move $m
157 > \rightarrow n$, and $\accMe$ is the probability of accepting the move
158 > $m \rightarrow n$.  Defining a probability vector,
159 > $\boldsymbol{\rho}$, such that
160 > \begin{equation}
161 > \boldsymbol{\rho} = \{\rho_1, \rho_2, \ldots \rho_m, \rho_n,
162 >        \ldots \rho_N \}
163 > \label{introEq:MCrhoVector}
164 > \end{equation}
165 > a transition matrix $\boldsymbol{\Pi}$ can be defined,
166 > whose elements are $\pi_{mn}$, for each given transition.  The
167 > limiting distribution of the Markov chain can then be found by
168 > applying the transition matrix an infinite number of times to the
169 > distribution vector.
170 > \begin{equation}
171 > \boldsymbol{\rho}_{\text{limit}} =
172 >        \lim_{N \rightarrow \infty} \boldsymbol{\rho}_{\text{initial}}
173 >        \boldsymbol{\Pi}^N
174 > \label{introEq:MCmarkovLimit}
175 > \end{equation}
176 > The limiting distribution of the chain is independent of the starting
177 > distribution, and successive applications of the transition matrix
178 > will only yield the limiting distribution again.
179 > \begin{equation}
180 > \boldsymbol{\rho}_{\text{limit}} = \boldsymbol{\rho}_{\text{initial}}
181 >        \boldsymbol{\Pi}
182 > \label{introEq:MCmarkovEquil}
183 > \end{equation}
184  
185 < symplectic methods
185 > \subsubsection{\label{introSec:metropolisMethod}The Metropolis Method}
186  
187 < Extended ensembles (NVT NPT)
187 > In the Metropolis method\cite{metropolis:1953}
188 > Eq.~\ref{introEq:MCmarkovEquil} is solved such that
189 > $\boldsymbol{\rho}_{\text{limit}}$ matches the Boltzman distribution
190 > of states.  The method accomplishes this by imposing the strong
191 > condition of microscopic reversibility on the equilibrium
192 > distribution.  Meaning, that at equilibrium the probability of going
193 > from $m$ to $n$ is the same as going from $n$ to $m$.
194 > \begin{equation}
195 > \rho_m\pi_{mn} = \rho_n\pi_{nm}
196 > \label{introEq:MCmicroReverse}
197 > \end{equation}
198 > Further, $\boldsymbol{\alpha}$ is chosen to be a symetric matrix in
199 > the Metropolis method.  Using Eq.~\ref{introEq:MCpi},
200 > Eq.~\ref{introEq:MCmicroReverse} becomes
201 > \begin{equation}
202 > \frac{\accMe(m \rightarrow n)}{\accMe(n \rightarrow m)} =
203 >        \frac{\rho_n}{\rho_m}
204 > \label{introEq:MCmicro2}
205 > \end{equation}
206 > For a Boltxman limiting distribution,
207 > \begin{equation}
208 > \frac{\rho_n}{\rho_m} = e^{-\beta[\mathcal{U}(n) - \mathcal{U}(m)]}
209 >        = e^{-\beta \Delta \mathcal{U}}
210 > \label{introEq:MCmicro3}
211 > \end{equation}
212 > This allows for the following set of acceptance rules be defined:
213 > \begin{equation}
214 > EQ Here
215 > \end{equation}
216  
217 < constrained dynamics
217 > Using the acceptance criteria from Eq.~\ref{fix} the Metropolis method
218 > proceeds as follows
219 > \begin{itemize}
220 > \item Generate an initial configuration $fix$ which has some finite probability in $fix$.
221 > \item Modify $fix$, to generate configuratioon $fix$.
222 > \item If configuration $n$ lowers the energy of the system, accept the move with unity ($fix$ becomes $fix$).  Otherwise accept with probability $fix$.
223 > \item Accumulate the average for the configurational observable of intereest.
224 > \item Repeat from step 2 until average converges.
225 > \end{itemize}
226 > One important note is that the average is accumulated whether the move
227 > is accepted or not, this ensures proper weighting of the average.
228 > Using Eq.~\ref{fix} it becomes clear that the accumulated averages are
229 > the ensemble averages, as this method ensures that the limiting
230 > distribution is the Boltzman distribution.
231  
232 + \subsection{\label{introSec:MD}Molecular Dynamics Simulations}
233 +
234 + The main simulation tool used in this research is Molecular Dynamics.
235 + Molecular Dynamics is when the equations of motion for a system are
236 + integrated in order to obtain information about both the positions and
237 + momentum of a system, allowing the calculation of not only
238 + configurational observables, but momenta dependent ones as well:
239 + diffusion constants, velocity auto correlations, folding/unfolding
240 + events, etc.  Due to the principle of ergodicity, Eq.~\ref{fix}, the
241 + average of these observables over the time period of the simulation
242 + are taken to be the ensemble averages for the system.
243 +
244 + The choice of when to use molecular dynamics over Monte Carlo
245 + techniques, is normally decided by the observables in which the
246 + researcher is interested.  If the observabvles depend on momenta in
247 + any fashion, then the only choice is molecular dynamics in some form.
248 + However, when the observable is dependent only on the configuration,
249 + then most of the time Monte Carlo techniques will be more efficent.
250 +
251 + The focus of research in the second half of this dissertation is
252 + centered around the dynamic properties of phospholipid bilayers,
253 + making molecular dynamics key in the simulation of those properties.
254 +
255 + \subsubsection{Molecular dynamics Algorithm}
256 +
257 + To illustrate how the molecular dynamics technique is applied, the
258 + following sections will describe the sequence involved in a
259 + simulation.  Sec.~\ref{fix} deals with the initialization of a
260 + simulation.  Sec.~\ref{fix} discusses issues involved with the
261 + calculation of the forces.  Sec.~\ref{fix} concludes the algorithm
262 + discussion with the integration of the equations of motion. \cite{fix}
263 +
264 + \subsubsection{initialization}
265 +
266 + When selecting the initial configuration for the simulation it is
267 + important to consider what dynamics one is hoping to observe.
268 + Ch.~\ref{fix} deals with the formation and equilibrium dynamics of
269 + phospholipid membranes.  Therefore in these simulations initial
270 + positions were selected that in some cases dispersed the lipids in
271 + water, and in other cases structured the lipids into preformed
272 + bilayers.  Important considerations at this stage of the simulation are:
273 + \begin{itemize}
274 + \item There are no major overlaps of molecular or atomic orbitals
275 + \item Velocities are chosen in such a way as to not gie the system a non=zero total momentum or angular momentum.
276 + \item It is also sometimes desireable to select the velocities to correctly sample the target temperature.
277 + \end{itemize}
278 +
279 + The first point is important due to the amount of potential energy
280 + generated by having two particles too close together.  If overlap
281 + occurs, the first evaluation of forces will return numbers so large as
282 + to render the numerical integration of teh motion meaningless.  The
283 + second consideration keeps the system from drifting or rotating as a
284 + whole.  This arises from the fact that most simulations are of systems
285 + in equilibrium in the absence of outside forces.  Therefore any net
286 + movement would be unphysical and an artifact of the simulation method
287 + used.  The final point addresses teh selection of the magnitude of the
288 + initial velocities.  For many simulations it is convienient to use
289 + this opportunity to scale the amount of kinetic energy to reflect the
290 + desired thermal distribution of the system.  However, it must be noted
291 + that most systems will require further velocity rescaling after the
292 + first few initial simulation steps due to either loss or gain of
293 + kinetic energy from energy stored in potential degrees of freedom.
294 +
295 + \subsubsection{Force Evaluation}
296 +
297 + The evaluation of forces is the most computationally expensive portion
298 + of a given molecular dynamics simulation.  This is due entirely to the
299 + evaluation of long range forces in a simulation, typically pair-wise.
300 + These forces are most commonly the Van der Waals force, and sometimes
301 + Coulombic forces as well.  For a pair-wise force, there are $fix$
302 + pairs to be evaluated, where $n$ is the number of particles in the
303 + system.  This leads to the calculations scaling as $fix$, making large
304 + simulations prohibitive in the absence of any computation saving
305 + techniques.
306 +
307 + Another consideration one must resolve, is that in a given simulation
308 + a disproportionate number of the particles will feel the effects of
309 + the surface. \cite{fix} For a cubic system of 1000 particles arranged
310 + in a $10x10x10$ cube, 488 particles will be exposed to the surface.
311 + Unless one is simulating an isolated particle group in a vacuum, the
312 + behavior of the system will be far from the desired bulk
313 + charecteristics.  To offset this, simulations employ the use of
314 + periodic boundary images. \cite{fix}
315 +
316 + The technique involves the use of an algorithm that replicates the
317 + simulation box on an infinite lattice in cartesian space.  Any given
318 + particle leaving the simulation box on one side will have an image of
319 + itself enter on the opposite side (see Fig.~\ref{fix}).
320 + \begin{equation}
321 + EQ Here
322 + \end{equation}
323 + In addition, this sets that any given particle pair has an image, real
324 + or periodic, within $fix$ of each other.  A discussion of the method
325 + used to calculate the periodic image can be found in Sec.\ref{fix}.
326 +
327 + Returning to the topic of the computational scale of the force
328 + evaluation, the use of periodic boundary conditions requires that a
329 + cutoff radius be employed.  Using a cutoff radius improves the
330 + efficiency of the force evaluation, as particles farther than a
331 + predetermined distance, $fix$, are not included in the
332 + calculation. \cite{fix} In a simultation with periodic images, $fix$
333 + has a maximum value of $fix$.  Fig.~\ref{fix} illustrates how using an
334 + $fix$ larger than this value, or in the extreme limit of no $fix$ at
335 + all, the corners of the simulation box are unequally weighted due to
336 + the lack of particle images in the $x$, $y$, or $z$ directions past a
337 + disance of $fix$.
338 +
339 + With the use of an $fix$, however, comes a discontinuity in the potential energy curve (Fig.~\ref{fix}).
340 +
341 +
342   \section{\label{introSec:chapterLayout}Chapter Layout}
343  
344   \subsection{\label{introSec:RSA}Random Sequential Adsorption}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines