| 1 |
|
|
| 2 |
|
|
| 3 |
< |
\chapter{\label{chapt:intro}Introduction and Theoretical Background} |
| 3 |
> |
\chapter{\label{chapt:intro}INTRODUCTION AND THEORETICAL BACKGROUND} |
| 4 |
|
|
| 5 |
|
|
| 6 |
|
The techniques used in the course of this research fall under the two |
| 7 |
|
main classes of molecular simulation: Molecular Dynamics and Monte |
| 8 |
< |
Carlo. Molecular Dynamic simulations integrate the equations of motion |
| 9 |
< |
for a given system of particles, allowing the researcher to gain |
| 10 |
< |
insight into the time dependent evolution of a system. Diffusion |
| 8 |
> |
Carlo. Molecular Dynamics simulations integrate the equations of |
| 9 |
> |
motion for a given system of particles, allowing the researcher to |
| 10 |
> |
gain insight into the time dependent evolution of a system. Diffusion |
| 11 |
|
phenomena are readily studied with this simulation technique, making |
| 12 |
|
Molecular Dynamics the main simulation technique used in this |
| 13 |
|
research. Other aspects of the research fall under the Monte Carlo |
| 14 |
|
class of simulations. In Monte Carlo, the configuration space |
| 15 |
< |
available to the collection of particles is sampled stochastically, |
| 16 |
< |
or randomly. Each configuration is chosen with a given probability |
| 17 |
< |
based on the Maxwell Boltzmann distribution. These types of simulations |
| 18 |
< |
are best used to probe properties of a system that are only dependent |
| 19 |
< |
only on the state of the system. Structural information about a system |
| 20 |
< |
is most readily obtained through these types of methods. |
| 15 |
> |
available to the collection of particles is sampled stochastically, or |
| 16 |
> |
randomly. Each configuration is chosen with a given probability based |
| 17 |
> |
on the Maxwell Boltzmann distribution. These types of simulations are |
| 18 |
> |
best used to probe properties of a system that are dependent only on |
| 19 |
> |
the state of the system. Structural information about a system is most |
| 20 |
> |
readily obtained through these types of methods. |
| 21 |
|
|
| 22 |
|
Although the two techniques employed seem dissimilar, they are both |
| 23 |
|
linked by the overarching principles of Statistical |
| 24 |
< |
Thermodynamics. Statistical Thermodynamics governs the behavior of |
| 24 |
> |
Mechanics. Statistical Meachanics governs the behavior of |
| 25 |
|
both classes of simulations and dictates what each method can and |
| 26 |
|
cannot do. When investigating a system, one most first analyze what |
| 27 |
|
thermodynamic properties of the system are being probed, then chose |
| 31 |
|
|
| 32 |
|
The following section serves as a brief introduction to some of the |
| 33 |
|
Statistical Mechanics concepts present in this dissertation. What |
| 34 |
< |
follows is a brief derivation of Boltzmann weighted statistics, and an |
| 34 |
> |
follows is a brief derivation of Boltzmann weighted statistics and an |
| 35 |
|
explanation of how one can use the information to calculate an |
| 36 |
|
observable for a system. This section then concludes with a brief |
| 37 |
|
discussion of the ergodic hypothesis and its relevance to this |
| 39 |
|
|
| 40 |
|
\subsection{\label{introSec:boltzman}Boltzmann weighted statistics} |
| 41 |
|
|
| 42 |
< |
Consider a system, $\gamma$, with some total energy,, $E_{\gamma}$. |
| 43 |
< |
Let $\Omega(E_{\gamma})$ represent the number of degenerate ways |
| 44 |
< |
$\boldsymbol{\Gamma}$, the collection of positions and conjugate |
| 45 |
< |
momenta of system $\gamma$, can be configured to give |
| 46 |
< |
$E_{\gamma}$. Further, if $\gamma$ is in contact with a bath system |
| 47 |
< |
where energy is exchanged between the two systems, $\Omega(E)$, where |
| 48 |
< |
$E$ is the total energy of both systems, can be represented as |
| 42 |
> |
Consider a system, $\gamma$, with total energy $E_{\gamma}$. Let |
| 43 |
> |
$\Omega(E_{\gamma})$ represent the number of degenerate ways |
| 44 |
> |
$\boldsymbol{\Gamma}\{r_1,r_2,\ldots r_n,p_1,p_2,\ldots p_n\}$, the |
| 45 |
> |
collection of positions and conjugate momenta of system $\gamma$, can |
| 46 |
> |
be configured to give $E_{\gamma}$. Further, if $\gamma$ is a subset |
| 47 |
> |
of a larger system, $\boldsymbol{\Lambda}\{E_1,E_2,\ldots |
| 48 |
> |
E_{\gamma},\ldots E_n\}$, the total degeneracy of the system can be |
| 49 |
> |
expressed as, |
| 50 |
|
\begin{equation} |
| 51 |
+ |
\Omega(\boldsymbol{\Lambda}) = \Omega(E_1) \times \Omega(E_2) \times \ldots |
| 52 |
+ |
\Omega(E_{\gamma}) \times \ldots \Omega(E_n) |
| 53 |
+ |
\label{introEq:SM0.1} |
| 54 |
+ |
\end{equation} |
| 55 |
+ |
This multiplicative combination of degeneracies is illustrated in |
| 56 |
+ |
Fig.~\ref{introFig:degenProd}. |
| 57 |
+ |
|
| 58 |
+ |
\begin{figure} |
| 59 |
+ |
\centering |
| 60 |
+ |
\includegraphics[width=\linewidth]{omegaFig.eps} |
| 61 |
+ |
\caption[An explanation of the combination of degeneracies]{Systems A and B both have three energy levels and two indistinguishable particles. When the total energy is 2, there are two ways for each system to disperse the energy. However, for system C, the superset of A and B, the total degeneracy is the product of the degeneracy of each system. In this case $\Omega(\text{C})$ is 4.} |
| 62 |
+ |
\label{introFig:degenProd} |
| 63 |
+ |
\end{figure} |
| 64 |
+ |
|
| 65 |
+ |
Next, consider the specific case of $\gamma$ in contact with a |
| 66 |
+ |
bath. Exchange of energy is allowed between the bath and the system, |
| 67 |
+ |
subject to the constraint that the total energy |
| 68 |
+ |
($E_{\text{bath}}+E_{\gamma}$) remain constant. $\Omega(E)$, where $E$ |
| 69 |
+ |
is the total energy of both systems, can be represented as |
| 70 |
+ |
\begin{equation} |
| 71 |
|
\Omega(E) = \Omega(E_{\gamma}) \times \Omega(E - E_{\gamma}) |
| 72 |
|
\label{introEq:SM1} |
| 73 |
|
\end{equation} |
| 78 |
|
\end{equation} |
| 79 |
|
|
| 80 |
|
The solution to Eq.~\ref{introEq:SM2} maximizes the number of |
| 81 |
< |
degenerative configurations in $E$. \cite{Frenkel1996} |
| 81 |
> |
degenerate configurations in $E$. \cite{Frenkel1996} |
| 82 |
|
This gives |
| 83 |
|
\begin{equation} |
| 84 |
|
\frac{\partial \ln \Omega(E)}{\partial E_{\gamma}} = 0 = |
| 197 |
|
|
| 198 |
|
\subsection{\label{introSec:ergodic}The Ergodic Hypothesis} |
| 199 |
|
|
| 200 |
< |
One last important consideration is that of ergodicity. Ergodicity is |
| 201 |
< |
the assumption that given an infinite amount of time, a system will |
| 202 |
< |
visit every available point in phase space.\cite{Frenkel1996} For most |
| 203 |
< |
systems, this is a valid assumption, except in cases where the system |
| 204 |
< |
may be trapped in a local feature (\emph{e.g.}~glasses). When valid, |
| 205 |
< |
ergodicity allows the unification of a time averaged observation and |
| 206 |
< |
an ensemble averaged one. If an observation is averaged over a |
| 207 |
< |
sufficiently long time, the system is assumed to visit all |
| 208 |
< |
appropriately available points in phase space, giving a properly |
| 209 |
< |
weighted statistical average. This allows the researcher freedom of |
| 210 |
< |
choice when deciding how best to measure a given observable. When an |
| 211 |
< |
ensemble averaged approach seems most logical, the Monte Carlo |
| 212 |
< |
techniques described in Sec.~\ref{introSec:monteCarlo} can be utilized. |
| 213 |
< |
Conversely, if a problem lends itself to a time averaging approach, |
| 214 |
< |
the Molecular Dynamics techniques in Sec.~\ref{introSec:MD} can be |
| 215 |
< |
employed. |
| 200 |
> |
In the case of a Molecular Dynamics simulation, rather than |
| 201 |
> |
calculating an ensemble average integral over phase space as in |
| 202 |
> |
Eq.~\ref{introEq:SM15}, it becomes easier to caclulate the time |
| 203 |
> |
average of an observable. Namely, |
| 204 |
> |
\begin{equation} |
| 205 |
> |
\langle A \rangle_t = \frac{1}{\tau} |
| 206 |
> |
\int_0^{\tau} A[\boldsymbol{\Gamma}(t)]\,dt |
| 207 |
> |
\label{introEq:SM16} |
| 208 |
> |
\end{equation} |
| 209 |
> |
Where the value of an observable is averaged over the length of time |
| 210 |
> |
that the simulation is run. This type of measurement mirrors the |
| 211 |
> |
experimental measurement of an observable. In an experiment, the |
| 212 |
> |
instrument analyzing the system must average its observation of the |
| 213 |
> |
finite time of the measurement. What is required then, is a principle |
| 214 |
> |
to relate the time average to the ensemble average. This is the |
| 215 |
> |
ergodic hypothesis. |
| 216 |
|
|
| 217 |
+ |
Ergodicity is the assumption that given an infinite amount of time, a |
| 218 |
+ |
system will visit every available point in phase |
| 219 |
+ |
space.\cite{Frenkel1996} For most systems, this is a valid assumption, |
| 220 |
+ |
except in cases where the system may be trapped in a local feature |
| 221 |
+ |
(\emph{e.g.}~glasses). When valid, ergodicity allows the unification |
| 222 |
+ |
of a time averaged observation and an ensemble averaged one. If an |
| 223 |
+ |
observation is averaged over a sufficiently long time, the system is |
| 224 |
+ |
assumed to visit all appropriately available points in phase space, |
| 225 |
+ |
giving a properly weighted statistical average. This allows the |
| 226 |
+ |
researcher freedom of choice when deciding how best to measure a given |
| 227 |
+ |
observable. When an ensemble averaged approach seems most logical, |
| 228 |
+ |
the Monte Carlo techniques described in Sec.~\ref{introSec:monteCarlo} |
| 229 |
+ |
can be utilized. Conversely, if a problem lends itself to a time |
| 230 |
+ |
averaging approach, the Molecular Dynamics techniques in |
| 231 |
+ |
Sec.~\ref{introSec:MD} can be employed. |
| 232 |
+ |
|
| 233 |
|
\section{\label{introSec:monteCarlo}Monte Carlo Simulations} |
| 234 |
|
|
| 235 |
|
The Monte Carlo method was developed by Metropolis and Ulam for their |
| 247 |
|
\end{equation} |
| 248 |
|
The equation can be recast as: |
| 249 |
|
\begin{equation} |
| 250 |
< |
I = (b-a)\langle f(x) \rangle |
| 250 |
> |
I = \int^b_a \frac{f(x)}{\rho(x)} \rho(x) dx |
| 251 |
> |
\label{introEq:Importance1} |
| 252 |
> |
\end{equation} |
| 253 |
> |
Where $\rho(x)$ is an arbitrary probability distribution in $x$. If |
| 254 |
> |
one conducts $\tau$ trials selecting a random number, $\zeta_\tau$, |
| 255 |
> |
from the distribution $\rho(x)$ on the interval $[a,b]$, then |
| 256 |
> |
Eq.~\ref{introEq:Importance1} becomes |
| 257 |
> |
\begin{equation} |
| 258 |
> |
I= \lim_{\tau \rightarrow \infty}\biggl \langle \frac{f(x)}{\rho(x)} \biggr \rangle_{\text{trials}[a,b]} |
| 259 |
> |
\label{introEq:Importance2} |
| 260 |
> |
\end{equation} |
| 261 |
> |
If $\rho(x)$ is uniformly distributed over the interval $[a,b]$, |
| 262 |
> |
\begin{equation} |
| 263 |
> |
\rho(x) = \frac{1}{b-a} |
| 264 |
> |
\label{introEq:importance2b} |
| 265 |
> |
\end{equation} |
| 266 |
> |
then the integral can be rewritten as |
| 267 |
> |
\begin{equation} |
| 268 |
> |
I = (b-a)\lim_{\tau \rightarrow \infty} |
| 269 |
> |
\langle f(x) \rangle_{\text{trials}[a,b]} |
| 270 |
|
\label{eq:MCex2} |
| 271 |
|
\end{equation} |
| 216 |
– |
Where $\langle f(x) \rangle$ is the unweighted average over the interval |
| 217 |
– |
$[a,b]$. The calculation of the integral could then be solved by |
| 218 |
– |
randomly choosing points along the interval $[a,b]$ and calculating |
| 219 |
– |
the value of $f(x)$ at each point. The accumulated average would then |
| 220 |
– |
approach $I$ in the limit where the number of trials is infinitely |
| 221 |
– |
large. |
| 272 |
|
|
| 273 |
|
However, in Statistical Mechanics, one is typically interested in |
| 274 |
|
integrals of the form: |
| 290 |
|
average. This is where importance sampling comes into |
| 291 |
|
play.\cite{allen87:csl} |
| 292 |
|
|
| 293 |
< |
Importance Sampling is a method where one selects a distribution from |
| 294 |
< |
which the random configurations are chosen in order to more |
| 295 |
< |
efficiently calculate the integral.\cite{Frenkel1996} Consider again |
| 296 |
< |
Eq.~\ref{eq:MCex1} rewritten to be: |
| 247 |
< |
\begin{equation} |
| 248 |
< |
I = \int^b_a \frac{f(x)}{\rho(x)} \rho(x) dx |
| 249 |
< |
\label{introEq:Importance1} |
| 250 |
< |
\end{equation} |
| 251 |
< |
Where $\rho(x)$ is an arbitrary probability distribution in $x$. If |
| 252 |
< |
one conducts $\tau$ trials selecting a random number, $\zeta_\tau$, |
| 253 |
< |
from the distribution $\rho(x)$ on the interval $[a,b]$, then |
| 254 |
< |
Eq.~\ref{introEq:Importance1} becomes |
| 255 |
< |
\begin{equation} |
| 256 |
< |
I= \biggl \langle \frac{f(x)}{\rho(x)} \biggr \rangle_{\text{trials}} |
| 257 |
< |
\label{introEq:Importance2} |
| 258 |
< |
\end{equation} |
| 259 |
< |
Looking at Eq.~\ref{eq:mcEnsAvg}, and realizing |
| 293 |
> |
Importance sampling is a method where the distribution, from which the |
| 294 |
> |
random configurations are chosen, is selected in such a way as to |
| 295 |
> |
efficiently sample the integral in question. Looking at |
| 296 |
> |
Eq.~\ref{eq:mcEnsAvg}, and realizing |
| 297 |
|
\begin {equation} |
| 298 |
|
\rho_{kT}(\mathbf{r}^N) = |
| 299 |
|
\frac{e^{-\beta V(\mathbf{r}^N)}} |
| 317 |
|
By selecting $\rho(\mathbf{r}^N)$ to be $\rho_{kT}(\mathbf{r}^N)$ |
| 318 |
|
Eq.~\ref{introEq:Importance4} becomes |
| 319 |
|
\begin{equation} |
| 320 |
< |
\langle A \rangle = \langle A(\mathbf{r}^N) \rangle_{\text{trials}} |
| 320 |
> |
\langle A \rangle = \langle A(\mathbf{r}^N) \rangle_{kT} |
| 321 |
|
\label{introEq:Importance5} |
| 322 |
|
\end{equation} |
| 323 |
|
The difficulty is selecting points $\mathbf{r}^N$ such that they are |
| 435 |
|
integrated in order to obtain information about both the positions and |
| 436 |
|
momentum of a system, allowing the calculation of not only |
| 437 |
|
configurational observables, but momenta dependent ones as well: |
| 438 |
< |
diffusion constants, velocity auto correlations, folding/unfolding |
| 439 |
< |
events, etc. Due to the principle of ergodicity, |
| 438 |
> |
diffusion constants, relaxation events, folding/unfolding |
| 439 |
> |
events, etc. With the time dependent information gained from a |
| 440 |
> |
Molecular Dynamics simulation, one can also calculate time correlation |
| 441 |
> |
functions of the form\cite{Hansen86} |
| 442 |
> |
\begin{equation} |
| 443 |
> |
\langle A(t)\,A(0)\rangle = \lim_{\tau\rightarrow\infty} \frac{1}{\tau} |
| 444 |
> |
\int_0^{\tau} A(t+t^{\prime})\,A(t^{\prime})\,dt^{\prime} |
| 445 |
> |
\label{introEq:timeCorr} |
| 446 |
> |
\end{equation} |
| 447 |
> |
These correlations can be used to measure fundamental time constants |
| 448 |
> |
of a system, such as diffusion constants from the velocity |
| 449 |
> |
autocorrelation or dipole relaxation times from the dipole |
| 450 |
> |
autocorrelation. Due to the principle of ergodicity, |
| 451 |
|
Sec.~\ref{introSec:ergodic}, the average of these observables over the |
| 452 |
|
time period of the simulation are taken to be the ensemble averages |
| 453 |
|
for the system. |
| 454 |
|
|
| 455 |
|
The choice of when to use molecular dynamics over Monte Carlo |
| 456 |
|
techniques, is normally decided by the observables in which the |
| 457 |
< |
researcher is interested. If the observables depend on momenta in |
| 457 |
> |
researcher is interested. If the observables depend on time in |
| 458 |
|
any fashion, then the only choice is molecular dynamics in some form. |
| 459 |
|
However, when the observable is dependent only on the configuration, |
| 460 |
< |
then most of the time Monte Carlo techniques will be more efficient. |
| 460 |
> |
then for most small systems, Monte Carlo techniques will be more efficient. |
| 461 |
|
|
| 462 |
|
The focus of research in the second half of this dissertation is |
| 463 |
|
centered around the dynamic properties of phospholipid bilayers, |
| 547 |
|
efficiency of the force evaluation, as particles farther than a |
| 548 |
|
predetermined distance, $r_{\text{cut}}$, are not included in the |
| 549 |
|
calculation.\cite{Frenkel1996} In a simulation with periodic images, |
| 550 |
< |
$r_{\text{cut}}$ has a maximum value of $\text{box}/2$. |
| 551 |
< |
Fig.~\ref{introFig:rMax} illustrates how when using an |
| 552 |
< |
$r_{\text{cut}}$ larger than this value, or in the extreme limit of no |
| 553 |
< |
$r_{\text{cut}}$ at all, the corners of the simulation box are |
| 554 |
< |
unequally weighted due to the lack of particle images in the $x$, $y$, |
| 555 |
< |
or $z$ directions past a distance of $\text{box} / 2$. |
| 550 |
> |
there are two methods to choose from, both with their own cutoff |
| 551 |
> |
limits. In the minimum image convention, $r_{\text{cut}}$ has a |
| 552 |
> |
maximum value of $\text{box}/2$. This is because each atom has only |
| 553 |
> |
one image that is seen by another atom, and further the image used is |
| 554 |
> |
the one that minimizes the distance between the two atoms. A system of |
| 555 |
> |
wrapped images about a central atom therefore has a maximum length |
| 556 |
> |
scale of box on a side (Fig.~\ref{introFig:rMaxMin}). The second |
| 557 |
> |
convention, multiple image convention, has a maximum $r_{\text{cut}}$ |
| 558 |
> |
of box. Here multiple images of each atom are replicated in the |
| 559 |
> |
periodic cells surrounding the central atom, this causes the atom to |
| 560 |
> |
see multiple copies of several atoms. If the cutoff radius is larger |
| 561 |
> |
than box, however, then the atom will see an image of itself, and |
| 562 |
> |
attempt to calculate an unphysical self-self force interaction |
| 563 |
> |
(Fig.~\ref{introFig:rMaxMult}). Due to the increased complexity and |
| 564 |
> |
commputaional ineffeciency of the multiple image method, the minimum |
| 565 |
> |
image method is the periodic method used throughout this research. |
| 566 |
|
|
| 567 |
|
\begin{figure} |
| 568 |
|
\centering |
| 569 |
|
\includegraphics[width=\linewidth]{rCutMaxFig.eps} |
| 570 |
< |
\caption[An explanation of $r_{\text{cut}}$]{The yellow atom has all other images wrapped to itself as the center. If $r_{\text{cut}}=\text{box}/2$, then the distribution is uniform (blue atoms). However, when $r_{\text{cut}}>\text{box}/2$ the corners are disproportionately weighted (green atoms) vs the axial directions (shaded regions).} |
| 571 |
< |
\label{introFig:rMax} |
| 570 |
> |
\caption[An explanation of minimum image convention]{The yellow atom has all other images wrapped to itself as the center. If $r_{\text{cut}}=\text{box}/2$, then the distribution is uniform (blue atoms). However, when $r_{\text{cut}}>\text{box}/2$ the corners are disproportionately weighted (green atoms) vs the axial directions (shaded regions).} |
| 571 |
> |
\label{introFig:rMaxMin} |
| 572 |
> |
\end{figure} |
| 573 |
> |
|
| 574 |
> |
\begin{figure} |
| 575 |
> |
\centering |
| 576 |
> |
\includegraphics[width=\linewidth]{rCutMaxMultFig.eps} |
| 577 |
> |
\caption[An explanation of multiple image convention]{The yellow atom is the central wrapping point. The blue atoms are the minimum images of the system about the central atom. The boxes with the green atoms are multiple images of the central box. If $r_{\text{cut}} \geq \{text{box}$ then the central atom sees multiple images of itself (red atom), creating a self-self force evaluation.} |
| 578 |
> |
\label{introFig:rMaxMult} |
| 579 |
|
\end{figure} |
| 580 |
|
|
| 581 |
|
With the use of an $r_{\text{cut}}$, however, comes a discontinuity in |
| 709 |
|
\begin{figure} |
| 710 |
|
\centering |
| 711 |
|
\includegraphics[width=\linewidth]{eulerRotFig.eps} |
| 712 |
< |
\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).} |
| 712 |
> |
\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$ axis (green rotation). Lastly is a final rotation of $\psi$ about the new $z$ axis (red rotation).} |
| 713 |
|
\label{introFig:eulerAngles} |
| 714 |
|
\end{figure} |
| 715 |
|
|
