| 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 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, 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. |
| 8 |
> |
Carlo. Molecular Dynamics simulations are carried out by integrating |
| 9 |
> |
the equations of motion for a given system of particles, allowing the |
| 10 |
> |
researcher to gain insight into the time dependent evolution of a |
| 11 |
> |
system. Transport phenomena are readily studied with this simulation |
| 12 |
> |
technique, making Molecular Dynamics the main simulation technique |
| 13 |
> |
used in this research. Other aspects of the research fall in the |
| 14 |
> |
Monte Carlo class of simulations. In Monte Carlo, the configuration |
| 15 |
> |
space available to the collection of particles is sampled |
| 16 |
> |
stochastically, or randomly. Each configuration is chosen with a given |
| 17 |
> |
probability based on the Boltzmann distribution. These types |
| 18 |
> |
of simulations are best used to probe properties of a system that are |
| 19 |
> |
dependent only on the state of the system. Structural information |
| 20 |
> |
about a system is most readily obtained through these types of |
| 21 |
> |
methods. |
| 22 |
|
|
| 23 |
|
Although the two techniques employed seem dissimilar, they are both |
| 24 |
|
linked by the overarching principles of Statistical |
| 25 |
< |
Mechanics. Statistical Meachanics governs the behavior of |
| 25 |
> |
Mechanics. Statistical Mechanics governs the behavior of |
| 26 |
|
both classes of simulations and dictates what each method can and |
| 27 |
< |
cannot do. When investigating a system, one most first analyze what |
| 28 |
< |
thermodynamic properties of the system are being probed, then chose |
| 27 |
> |
cannot do. When investigating a system, one must first analyze what |
| 28 |
> |
thermodynamic properties of the system are being probed, then choose |
| 29 |
|
which method best suits that objective. |
| 30 |
|
|
| 31 |
|
\section{\label{introSec:statThermo}Statistical Mechanics} |
| 32 |
|
|
| 33 |
|
The following section serves as a brief introduction to some of the |
| 34 |
< |
Statistical Mechanics concepts present in this dissertation. What |
| 34 |
> |
Statistical Mechanics concepts presented in this dissertation. What |
| 35 |
|
follows is a brief derivation of Boltzmann weighted statistics and an |
| 36 |
|
explanation of how one can use the information to calculate an |
| 37 |
|
observable for a system. This section then concludes with a brief |
| 125 |
|
\beta( E_{\gamma} ) = \beta( E_{\text{bath}} ) |
| 126 |
|
\label{introEq:SM8} |
| 127 |
|
\end{equation} |
| 128 |
< |
Showing that the partitioning of energy between the two systems is |
| 129 |
< |
actually a process of thermal equilibration.\cite{Frenkel1996} |
| 128 |
> |
Eq.~\ref{introEq:SM8} shows that the partitioning of energy between |
| 129 |
> |
the two systems is actually a process of thermal |
| 130 |
> |
equilibration.\cite{Frenkel1996} |
| 131 |
|
|
| 132 |
< |
An application of these results is to formulate the form of an |
| 133 |
< |
expectation value of an observable, $A$, in the canonical ensemble. In |
| 134 |
< |
the canonical ensemble, the number of particles, $N$, the volume, $V$, |
| 135 |
< |
and the temperature, $T$, are all held constant while the energy, $E$, |
| 136 |
< |
is allowed to fluctuate. Returning to the previous example, the bath |
| 132 |
> |
An application of these results is to formulate the expectation value |
| 133 |
> |
of an observable, $A$, in the canonical ensemble. In the canonical |
| 134 |
> |
ensemble, the number of particles, $N$, the volume, $V$, and the |
| 135 |
> |
temperature, $T$, are all held constant while the energy, $E$, is |
| 136 |
> |
allowed to fluctuate. Returning to the previous example, the bath |
| 137 |
|
system is now an infinitely large thermal bath, whose exchange of |
| 138 |
|
energy with the system $\gamma$ holds the temperature constant. The |
| 139 |
< |
partitioning of energy in the bath system is then related to the total |
| 140 |
< |
energy of both systems and the fluctuations in $E_{\gamma}$: |
| 139 |
> |
partitioning of energy between the bath and the system is then related |
| 140 |
> |
to the total energy of both systems and the fluctuations in |
| 141 |
> |
$E_{\gamma}$: |
| 142 |
|
\begin{equation} |
| 143 |
|
\Omega( E_{\gamma} ) = \Omega( E - E_{\gamma} ) |
| 144 |
|
\label{introEq:SM9} |
| 151 |
|
\label{introEq:SM10} |
| 152 |
|
\end{equation} |
| 153 |
|
Where $\int\limits_{\boldsymbol{\Gamma}} d\boldsymbol{\Gamma}$ denotes |
| 154 |
< |
an integration over all accessible phase space, $P_{\gamma}$ is the |
| 155 |
< |
probability of being in a given phase state and |
| 156 |
< |
$A(\boldsymbol{\Gamma})$ is some observable that is a function of the |
| 154 |
> |
an integration over all accessible points in phase space, $P_{\gamma}$ |
| 155 |
> |
is the probability of being in a given phase state and |
| 156 |
> |
$A(\boldsymbol{\Gamma})$ is an observable that is a function of the |
| 157 |
|
phase state. |
| 158 |
|
|
| 159 |
|
Because entropy seeks to maximize the number of degenerate states at a |
| 209 |
|
\int_0^{\tau} A[\boldsymbol{\Gamma}(t)]\,dt |
| 210 |
|
\label{introEq:SM16} |
| 211 |
|
\end{equation} |
| 212 |
< |
Where the value of an observable is averaged over the length of time |
| 213 |
< |
that the simulation is run. This type of measurement mirrors the |
| 214 |
< |
experimental measurement of an observable. In an experiment, the |
| 212 |
> |
Where the value of an observable is averaged over the length of time, |
| 213 |
> |
$\tau$, that the simulation is run. This type of measurement mirrors |
| 214 |
> |
the experimental measurement of an observable. In an experiment, the |
| 215 |
|
instrument analyzing the system must average its observation over the |
| 216 |
|
finite time of the measurement. What is required then, is a principle |
| 217 |
|
to relate the time average to the ensemble average. This is the |
| 224 |
|
(\emph{e.g.}~glasses). When valid, ergodicity allows the unification |
| 225 |
|
of a time averaged observation and an ensemble averaged one. If an |
| 226 |
|
observation is averaged over a sufficiently long time, the system is |
| 227 |
< |
assumed to visit all appropriately available points in phase space, |
| 227 |
> |
assumed to visit all energetically accessible points in phase space, |
| 228 |
|
giving a properly weighted statistical average. This allows the |
| 229 |
|
researcher freedom of choice when deciding how best to measure a given |
| 230 |
|
observable. When an ensemble averaged approach seems most logical, |
| 237 |
|
|
| 238 |
|
The Monte Carlo method was developed by Metropolis and Ulam for their |
| 239 |
|
work in fissionable material.\cite{metropolis:1949} The method is so |
| 240 |
< |
named, because it heavily uses random numbers in its |
| 241 |
< |
solution.\cite{allen87:csl} The Monte Carlo method allows for the |
| 242 |
< |
solution of integrals through the stochastic sampling of the values |
| 243 |
< |
within the integral. In the simplest case, the evaluation of an |
| 244 |
< |
integral would follow a brute force method of |
| 245 |
< |
sampling.\cite{Frenkel1996} Consider the following single dimensional |
| 243 |
< |
integral: |
| 240 |
> |
named, because it uses random numbers extensively.\cite{allen87:csl} |
| 241 |
> |
The Monte Carlo method allows for the solution of integrals through |
| 242 |
> |
the stochastic sampling of the values within the integral. In the |
| 243 |
> |
simplest case, the evaluation of an integral would follow a brute |
| 244 |
> |
force method of sampling.\cite{Frenkel1996} Consider the following |
| 245 |
> |
single dimensional integral: |
| 246 |
|
\begin{equation} |
| 247 |
< |
I = f(x)dx |
| 247 |
> |
I = \int_a^b f(x)dx |
| 248 |
|
\label{eq:MCex1} |
| 249 |
|
\end{equation} |
| 250 |
|
The equation can be recast as: |
| 373 |
|
distribution, and successive applications of the transition matrix |
| 374 |
|
will only yield the limiting distribution again. |
| 375 |
|
\begin{equation} |
| 376 |
< |
\boldsymbol{\rho}_{\text{limit}} = \boldsymbol{\rho}_{\text{initial}} |
| 376 |
> |
\boldsymbol{\rho}_{\text{limit}} = \boldsymbol{\rho}_{\text{limit}} |
| 377 |
|
\boldsymbol{\Pi} |
| 378 |
|
\label{introEq:MCmarkovEquil} |
| 379 |
|
\end{equation} |
| 385 |
|
$\boldsymbol{\rho}_{\text{limit}}$ matches the Boltzmann distribution |
| 386 |
|
of states. The method accomplishes this by imposing the strong |
| 387 |
|
condition of microscopic reversibility on the equilibrium |
| 388 |
< |
distribution. Meaning, that at equilibrium the probability of going |
| 388 |
> |
distribution. This means that at equilibrium, the probability of going |
| 389 |
|
from $m$ to $n$ is the same as going from $n$ to $m$. |
| 390 |
|
\begin{equation} |
| 391 |
|
\rho_m\pi_{mn} = \rho_n\pi_{nm} |
| 420 |
|
\begin{enumerate} |
| 421 |
|
\item Generate an initial configuration $\mathbf{r}^N$ which has some finite probability in $\rho_{kT}$. |
| 422 |
|
\item Modify $\mathbf{r}^N$, to generate configuration $\mathbf{r^{\prime}}^N$. |
| 423 |
< |
\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}}$. |
| 423 |
> |
\item If the new configuration lowers the energy of the system, accept the move ($\mathbf{r}^N$ becomes $\mathbf{r^{\prime}}^N$). Otherwise accept with probability $e^{-\beta \Delta \mathcal{U}}$. |
| 424 |
|
\item Accumulate the average for the configurational observable of interest. |
| 425 |
|
\item Repeat from step 2 until the average converges. |
| 426 |
|
\end{enumerate} |
| 433 |
|
\section{\label{introSec:MD}Molecular Dynamics Simulations} |
| 434 |
|
|
| 435 |
|
The main simulation tool used in this research is Molecular Dynamics. |
| 436 |
< |
Molecular Dynamics is when the equations of motion for a system are |
| 437 |
< |
integrated in order to obtain information about both the positions and |
| 438 |
< |
momentum of a system, allowing the calculation of not only |
| 439 |
< |
configurational observables, but momenta dependent ones as well: |
| 440 |
< |
diffusion constants, relaxation events, folding/unfolding |
| 441 |
< |
events, etc. With the time dependent information gained from a |
| 442 |
< |
Molecular Dynamics simulation, one can also calculate time correlation |
| 443 |
< |
functions of the form\cite{Hansen86} |
| 436 |
> |
Molecular Dynamics is the numerical integration of the equations of |
| 437 |
> |
motion for a system in order to obtain information about the dynamic |
| 438 |
> |
changes in the positions and momentum of the constituent particles. |
| 439 |
> |
This allows the calculation of not only configurational observables, |
| 440 |
> |
but momentum dependent ones as well: diffusion constants, relaxation |
| 441 |
> |
events, folding/unfolding events, etc. With the time dependent |
| 442 |
> |
information gained from a Molecular Dynamics simulation, one can also |
| 443 |
> |
calculate time correlation functions of the form\cite{Hansen86} |
| 444 |
|
\begin{equation} |
| 445 |
|
\langle A(t)\,A(0)\rangle = \lim_{\tau\rightarrow\infty} \frac{1}{\tau} |
| 446 |
|
\int_0^{\tau} A(t+t^{\prime})\,A(t^{\prime})\,dt^{\prime} |
| 450 |
|
of a system, such as diffusion constants from the velocity |
| 451 |
|
autocorrelation or dipole relaxation times from the dipole |
| 452 |
|
autocorrelation. Due to the principle of ergodicity, |
| 453 |
< |
Sec.~\ref{introSec:ergodic}, the average of these observables over the |
| 454 |
< |
time period of the simulation are taken to be the ensemble averages |
| 455 |
< |
for the system. |
| 453 |
> |
Sec.~\ref{introSec:ergodic}, the average of static observables over |
| 454 |
> |
the length of the simulation are taken to be the ensemble averages for |
| 455 |
> |
the system. |
| 456 |
|
|
| 457 |
|
The choice of when to use molecular dynamics over Monte Carlo |
| 458 |
|
techniques, is normally decided by the observables in which the |
| 459 |
< |
researcher is interested. If the observables depend on time in |
| 460 |
< |
any fashion, then the only choice is molecular dynamics in some form. |
| 459 |
> |
researcher is interested. If the observables depend on time in any |
| 460 |
> |
fashion, then the only choice is molecular dynamics in some form. |
| 461 |
|
However, when the observable is dependent only on the configuration, |
| 462 |
< |
then for most small systems, Monte Carlo techniques will be more efficient. |
| 462 |
> |
then for most small systems, Monte Carlo techniques will be more |
| 463 |
> |
efficient. The focus of research in the second half of this |
| 464 |
> |
dissertation is centered on the dynamic properties of phospholipid |
| 465 |
> |
bilayers, making molecular dynamics key in the simulation of those |
| 466 |
> |
properties. |
| 467 |
|
|
| 468 |
< |
The focus of research in the second half of this dissertation is |
| 463 |
< |
centered around the dynamic properties of phospholipid bilayers, |
| 464 |
< |
making molecular dynamics key in the simulation of those properties. |
| 468 |
> |
\subsection{\label{introSec:mdAlgorithm}Molecular Dynamics Algorithms} |
| 469 |
|
|
| 466 |
– |
\subsection{\label{introSec:mdAlgorithm}The Molecular Dynamics Algorithm} |
| 467 |
– |
|
| 470 |
|
To illustrate how the molecular dynamics technique is applied, the |
| 471 |
|
following sections will describe the sequence involved in a |
| 472 |
|
simulation. Sec.~\ref{introSec:mdInit} deals with the initialization |
| 482 |
|
Ch.~\ref{chapt:lipid} deals with the formation and equilibrium dynamics of |
| 483 |
|
phospholipid membranes. Therefore in these simulations initial |
| 484 |
|
positions were selected that in some cases dispersed the lipids in |
| 485 |
< |
water, and in other cases structured the lipids into performed |
| 485 |
> |
water, and in other cases structured the lipids into preformed |
| 486 |
|
bilayers. Important considerations at this stage of the simulation are: |
| 487 |
|
\begin{itemize} |
| 488 |
< |
\item There are no major overlaps of molecular or atomic orbitals |
| 489 |
< |
\item Velocities are chosen in such a way as to not give the system a non=zero total momentum or angular momentum. |
| 490 |
< |
\item It is also sometimes desirable to select the velocities to correctly sample the target temperature. |
| 488 |
> |
\item There are no overlapping molecules or atoms. |
| 489 |
> |
\item Velocities are chosen in such a way as to give the system zero total momentum and angular momentum. |
| 490 |
> |
\item It is also sometimes desirable to select the velocities to correctly sample the ensemble at a particular target temperature. |
| 491 |
|
\end{itemize} |
| 492 |
|
|
| 493 |
< |
The first point is important due to the amount of potential energy |
| 494 |
< |
generated by having two particles too close together. If overlap |
| 495 |
< |
occurs, the first evaluation of forces will return numbers so large as |
| 496 |
< |
to render the numerical integration of the motion meaningless. The |
| 497 |
< |
second consideration keeps the system from drifting or rotating as a |
| 498 |
< |
whole. This arises from the fact that most simulations are of systems |
| 499 |
< |
in equilibrium in the absence of outside forces. Therefore any net |
| 493 |
> |
The first point is important due to the forces generated when two |
| 494 |
> |
particles are too close together. If overlap occurs, the first |
| 495 |
> |
evaluation of forces will return numbers so large as to render the |
| 496 |
> |
numerical integration of the equations motion meaningless. The second |
| 497 |
> |
consideration keeps the system from drifting or rotating as a whole. |
| 498 |
> |
This arises from the fact that most simulations are of systems in |
| 499 |
> |
equilibrium in the absence of outside forces. Therefore any net |
| 500 |
|
movement would be unphysical and an artifact of the simulation method |
| 501 |
|
used. The final point addresses the selection of the magnitude of the |
| 502 |
< |
initial velocities. For many simulations it is convenient to use |
| 503 |
< |
this opportunity to scale the amount of kinetic energy to reflect the |
| 502 |
> |
initial velocities. For many simulations it is convenient to use this |
| 503 |
> |
opportunity to scale the amount of kinetic energy to reflect the |
| 504 |
|
desired thermal distribution of the system. However, it must be noted |
| 505 |
< |
that most systems will require further velocity rescaling after the |
| 506 |
< |
first few initial simulation steps due to either loss or gain of |
| 507 |
< |
kinetic energy from energy stored in potential degrees of freedom. |
| 505 |
> |
that due to equipartition, most systems will require further velocity |
| 506 |
> |
rescaling after the first few initial simulation steps due to either |
| 507 |
> |
loss or gain of kinetic energy from energy stored as potential energy. |
| 508 |
|
|
| 509 |
|
\subsection{\label{introSec:mdForce}Force Evaluation} |
| 510 |
|
|
| 511 |
|
The evaluation of forces is the most computationally expensive portion |
| 512 |
< |
of a given molecular dynamics simulation. This is due entirely to the |
| 513 |
< |
evaluation of long range forces in a simulation, typically pair-wise. |
| 514 |
< |
These forces are most commonly the Van der Waals force, and sometimes |
| 515 |
< |
Coulombic forces as well. For a pair-wise force, there are $N(N-1)/ 2$ |
| 516 |
< |
pairs to be evaluated, where $N$ is the number of particles in the |
| 517 |
< |
system. This leads to the calculations scaling as $N^2$, making large |
| 518 |
< |
simulations prohibitive in the absence of any computation saving |
| 519 |
< |
techniques. |
| 512 |
> |
of any molecular dynamics simulation. This is due entirely to the |
| 513 |
> |
evaluation of long range forces in a simulation, which are typically |
| 514 |
> |
pair potentials. These forces are most commonly the van der Waals |
| 515 |
> |
force, and sometimes Coulombic forces as well. For a pair-wise |
| 516 |
> |
interaction, there are $N(N-1)/ 2$ pairs to be evaluated, where $N$ is |
| 517 |
> |
the number of particles in the system. This leads to calculations which |
| 518 |
> |
scale as $N^2$, making large simulations prohibitive in the absence |
| 519 |
> |
of any computation saving techniques. |
| 520 |
|
|
| 521 |
< |
Another consideration one must resolve, is that in a given simulation |
| 521 |
> |
Another consideration one must resolve, is that in a given simulation, |
| 522 |
|
a disproportionate number of the particles will feel the effects of |
| 523 |
|
the surface.\cite{allen87:csl} For a cubic system of 1000 particles |
| 524 |
|
arranged in a $10 \times 10 \times 10$ cube, 488 particles will be |
| 525 |
< |
exposed to the surface. Unless one is simulating an isolated particle |
| 526 |
< |
group in a vacuum, the behavior of the system will be far from the |
| 527 |
< |
desired bulk characteristics. To offset this, simulations employ the |
| 528 |
< |
use of periodic boundary images.\cite{born:1912} |
| 525 |
> |
exposed to the surface. Unless one is simulating an isolated cluster |
| 526 |
> |
in a vacuum, the behavior of the system will be far from the desired |
| 527 |
> |
bulk characteristics. To offset this, simulations employ the use of |
| 528 |
> |
periodic images.\cite{born:1912} |
| 529 |
|
|
| 530 |
|
The technique involves the use of an algorithm that replicates the |
| 531 |
< |
simulation box on an infinite lattice in Cartesian space. Any given |
| 531 |
> |
simulation box on an infinite lattice in Cartesian space. Any |
| 532 |
|
particle leaving the simulation box on one side will have an image of |
| 533 |
< |
itself enter on the opposite side (see Fig.~\ref{introFig:pbc}). In |
| 534 |
< |
addition, this sets that any two particles have an image, real or |
| 535 |
< |
periodic, within $\text{box}/2$ of each other. A discussion of the |
| 536 |
< |
method used to calculate the periodic image can be found in |
| 537 |
< |
Sec.\ref{oopseSec:pbc}. |
| 533 |
> |
itself enter on the opposite side (see |
| 534 |
> |
Fig.~\ref{introFig:pbc}). Therefore, a pair of particles have an |
| 535 |
> |
image of each other, real or periodic, within $\text{box}/2$. A |
| 536 |
> |
discussion of the method used to calculate the periodic image can be |
| 537 |
> |
found in Sec.\ref{oopseSec:pbc}. |
| 538 |
|
|
| 539 |
|
\begin{figure} |
| 540 |
|
\centering |
| 543 |
|
\label{introFig:pbc} |
| 544 |
|
\end{figure} |
| 545 |
|
|
| 546 |
< |
Returning to the topic of the computational scale of the force |
| 546 |
> |
Returning to the topic of the computational scaling of the force |
| 547 |
|
evaluation, the use of periodic boundary conditions requires that a |
| 548 |
|
cutoff radius be employed. Using a cutoff radius improves the |
| 549 |
|
efficiency of the force evaluation, as particles farther than a |
| 552 |
|
there are two methods to choose from, both with their own cutoff |
| 553 |
|
limits. In the minimum image convention, $r_{\text{cut}}$ has a |
| 554 |
|
maximum value of $\text{box}/2$. This is because each atom has only |
| 555 |
< |
one image that is seen by another atom, and further the image used is |
| 556 |
< |
the one that minimizes the distance between the two atoms. A system of |
| 557 |
< |
wrapped images about a central atom therefore has a maximum length |
| 558 |
< |
scale of box on a side (Fig.~\ref{introFig:rMaxMin}). The second |
| 559 |
< |
convention, multiple image convention, has a maximum $r_{\text{cut}}$ |
| 560 |
< |
of box. Here multiple images of each atom are replicated in the |
| 555 |
> |
one image that is seen by other atoms. The image used is the one that |
| 556 |
> |
minimizes the distance between the two atoms. A system of wrapped |
| 557 |
> |
images about a central atom therefore has a maximum length scale of |
| 558 |
> |
box on a side (Fig.~\ref{introFig:rMaxMin}). The second convention, |
| 559 |
> |
multiple image convention, has a maximum $r_{\text{cut}}$ of the box |
| 560 |
> |
length itself. Here multiple images of each atom are replicated in the |
| 561 |
|
periodic cells surrounding the central atom, this causes the atom to |
| 562 |
|
see multiple copies of several atoms. If the cutoff radius is larger |
| 563 |
|
than box, however, then the atom will see an image of itself, and |
| 564 |
|
attempt to calculate an unphysical self-self force interaction |
| 565 |
|
(Fig.~\ref{introFig:rMaxMult}). Due to the increased complexity and |
| 566 |
|
commputaional ineffeciency of the multiple image method, the minimum |
| 567 |
< |
image method is the periodic method used throughout this research. |
| 567 |
> |
image method is the used throughout this research. |
| 568 |
|
|
| 569 |
|
\begin{figure} |
| 570 |
|
\centering |
| 576 |
|
\begin{figure} |
| 577 |
|
\centering |
| 578 |
|
\includegraphics[width=\linewidth]{rCutMaxMultFig.eps} |
| 579 |
< |
\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.} |
| 579 |
> |
\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.} |
| 580 |
|
\label{introFig:rMaxMult} |
| 581 |
|
\end{figure} |
| 582 |
|
|
| 583 |
< |
With the use of an $r_{\text{cut}}$, however, comes a discontinuity in |
| 583 |
> |
With the use of a cutoff radius, however, comes a discontinuity in |
| 584 |
|
the potential energy curve (Fig.~\ref{introFig:shiftPot}). To fix this |
| 585 |
|
discontinuity, one calculates the potential energy at the |
| 586 |
|
$r_{\text{cut}}$, and adds that value to the potential, causing |
| 596 |
|
\end{figure} |
| 597 |
|
|
| 598 |
|
The second main simplification used in this research is the Verlet |
| 599 |
< |
neighbor list. \cite{allen87:csl} In the Verlet method, one generates |
| 599 |
> |
neighbor list.\cite{allen87:csl} In the Verlet method, one generates |
| 600 |
|
a list of all neighbor atoms, $j$, surrounding atom $i$ within some |
| 601 |
|
cutoff $r_{\text{list}}$, where $r_{\text{list}}>r_{\text{cut}}$. |
| 602 |
|
This list is created the first time forces are evaluated, then on |
| 603 |
|
subsequent force evaluations, pair calculations are only calculated |
| 604 |
< |
from the neighbor lists. The lists are updated if any given particle |
| 604 |
> |
from the neighbor lists. The lists are updated if any particle |
| 605 |
|
in the system moves farther than $r_{\text{list}}-r_{\text{cut}}$, |
| 606 |
< |
giving rise to the possibility that a particle has left or joined a |
| 606 |
> |
which indeicates the possibility that a particle has left or joined the |
| 607 |
|
neighbor list. |
| 608 |
|
|
| 609 |
< |
\subsection{\label{introSec:mdIntegrate} Integration of the equations of motion} |
| 609 |
> |
\subsection{\label{introSec:mdIntegrate} Integration of the Equations of Motion} |
| 610 |
|
|
| 611 |
|
A starting point for the discussion of molecular dynamics integrators |
| 612 |
|
is the Verlet algorithm.\cite{Frenkel1996} It begins with a Taylor |
| 613 |
|
expansion of position in time: |
| 614 |
|
\begin{equation} |
| 615 |
|
q(t+\Delta t)= q(t) + v(t)\Delta t + \frac{F(t)}{2m}\Delta t^2 + |
| 616 |
< |
\frac{\Delta t^3}{3!}\frac{\partial q(t)}{\partial t} + |
| 616 |
> |
\frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} + |
| 617 |
|
\mathcal{O}(\Delta t^4) |
| 618 |
|
\label{introEq:verletForward} |
| 619 |
|
\end{equation} |
| 620 |
|
As well as, |
| 621 |
|
\begin{equation} |
| 622 |
|
q(t-\Delta t)= q(t) - v(t)\Delta t + \frac{F(t)}{2m}\Delta t^2 - |
| 623 |
< |
\frac{\Delta t^3}{3!}\frac{\partial q(t)}{\partial t} + |
| 623 |
> |
\frac{\Delta t^3}{3!}\frac{\partial^3 q(t)}{\partial t^3} + |
| 624 |
|
\mathcal{O}(\Delta t^4) |
| 625 |
|
\label{introEq:verletBack} |
| 626 |
|
\end{equation} |
| 643 |
|
order of $\Delta t^4$. |
| 644 |
|
|
| 645 |
|
In practice, however, the simulations in this research were integrated |
| 646 |
< |
with a velocity reformulation of the Verlet method.\cite{allen87:csl} |
| 646 |
> |
with the velocity reformulation of the Verlet method.\cite{allen87:csl} |
| 647 |
|
\begin{align} |
| 648 |
|
q(t+\Delta t) &= q(t) + v(t)\Delta t + \frac{F(t)}{2m}\Delta t^2 % |
| 649 |
|
\label{introEq:MDvelVerletPos} \\% |
| 654 |
|
The original Verlet algorithm can be regained by substituting the |
| 655 |
|
velocity back into Eq.~\ref{introEq:MDvelVerletPos}. The Verlet |
| 656 |
|
formulations are chosen in this research because the algorithms have |
| 657 |
< |
very little long term drift in energy conservation. Energy |
| 658 |
< |
conservation in a molecular dynamics simulation is of extreme |
| 659 |
< |
importance, as it is a measure of how closely one is following the |
| 660 |
< |
``true'' trajectory with the finite integration scheme. An exact |
| 661 |
< |
solution to the integration will conserve area in phase space, as well |
| 657 |
> |
very little long time drift in energy. Energy conservation in a |
| 658 |
> |
molecular dynamics simulation is of extreme importance, as it is a |
| 659 |
> |
measure of how closely one is following the ``true'' trajectory with |
| 660 |
> |
the finite integration scheme. An exact solution to the integration |
| 661 |
> |
will conserve area in phase space (i.e.~will be symplectic), as well |
| 662 |
|
as be reversible in time, that is, the trajectory integrated forward |
| 663 |
|
or backwards will exactly match itself. Having a finite algorithm |
| 664 |
|
that both conserves area in phase space and is time reversible, |
| 665 |
< |
therefore increases, but does not guarantee the ``correctness'' or the |
| 666 |
< |
integrated trajectory. |
| 665 |
> |
therefore increases the reliability, but does not guarantee the |
| 666 |
> |
``correctness'' of the integrated trajectory. |
| 667 |
|
|
| 668 |
|
It can be shown,\cite{Frenkel1996} that although the Verlet algorithm |
| 669 |
|
does not rigorously preserve the actual Hamiltonian, it does preserve |
| 679 |
|
\subsection{\label{introSec:MDfurther}Further Considerations} |
| 680 |
|
|
| 681 |
|
In the simulations presented in this research, a few additional |
| 682 |
< |
parameters are needed to describe the motions. The simulations |
| 683 |
< |
involving water and phospholipids in Ch.~\ref{chapt:lipid} are |
| 682 |
> |
parameters are needed to describe the motions. In the simulations |
| 683 |
> |
involving water and phospholipids in Ch.~\ref{chapt:lipid}, we are |
| 684 |
|
required to integrate the equations of motions for dipoles on atoms. |
| 685 |
|
This involves an additional three parameters be specified for each |
| 686 |
|
dipole atom: $\phi$, $\theta$, and $\psi$. These three angles are |
| 688 |
|
$z$-axis, and $\theta$ is a rotation about the new $x$-axis, and |
| 689 |
|
$\psi$ is a final rotation about the new $z$-axis (see |
| 690 |
|
Fig.~\ref{introFig:eulerAngles}). This sequence of rotations can be |
| 691 |
< |
accumulated into a single $3 \times 3$ matrix, $\mathbf{A}$, |
| 691 |
> |
accumulated into a single $3 \times 3$ matrix, $\mathsf{A}$, |
| 692 |
|
defined as follows: |
| 693 |
|
\begin{equation} |
| 694 |
< |
\mathbf{A} = |
| 694 |
> |
\mathsf{A} = |
| 695 |
|
\begin{bmatrix} |
| 696 |
|
\cos\phi\cos\psi-\sin\phi\cos\theta\sin\psi &% |
| 697 |
|
\sin\phi\cos\psi+\cos\phi\cos\theta\sin\psi &% |
| 730 |
|
\omega^s_y \frac{\cos\phi}{\sin\theta} |
| 731 |
|
\label{introEq:MDeulerPsi} |
| 732 |
|
\end{align} |
| 733 |
< |
Where $\omega^s_i$ is the angular velocity in the lab space frame |
| 734 |
< |
along Cartesian coordinate $i$. However, a difficulty arises when |
| 733 |
> |
Where $\omega^s_{\alpha}$ is the angular velocity in the lab space frame |
| 734 |
> |
along Cartesian coordinate $\alpha$. However, a difficulty arises when |
| 735 |
|
attempting to integrate Eq.~\ref{introEq:MDeulerPhi} and |
| 736 |
|
Eq.~\ref{introEq:MDeulerPsi}. The $\frac{1}{\sin \theta}$ present in |
| 737 |
|
both equations means there is a non-physical instability present when |
| 738 |
|
$\theta$ is 0 or $\pi$. To correct for this, the simulations integrate |
| 739 |
< |
the rotation matrix, $\mathbf{A}$, directly, thus avoiding the |
| 739 |
> |
the rotation matrix, $\mathsf{A}$, directly, thus avoiding the |
| 740 |
|
instability. This method was proposed by Dullweber |
| 741 |
|
\emph{et. al.}\cite{Dullweber1997}, and is presented in |
| 742 |
|
Sec.~\ref{introSec:MDsymplecticRot}. |
| 745 |
|
|
| 746 |
|
Before discussing the integration of the rotation matrix, it is |
| 747 |
|
necessary to understand the construction of a ``good'' integration |
| 748 |
< |
scheme. It has been previously |
| 749 |
< |
discussed(Sec.~\ref{introSec:mdIntegrate}) how it is desirable for an |
| 750 |
< |
integrator to be symplectic, or time reversible. The following is an |
| 748 |
> |
scheme. It has been previously discussed |
| 749 |
> |
(Sec.~\ref{introSec:mdIntegrate}) how it is desirable for an |
| 750 |
> |
integrator to be symplectic and time reversible. The following is an |
| 751 |
|
outline of the Trotter factorization of the Liouville Propagator as a |
| 752 |
< |
scheme for generating symplectic integrators.\cite{Tuckerman92} |
| 752 |
> |
scheme for generating symplectic, time-reversible |
| 753 |
> |
integrators.\cite{Tuckerman92} |
| 754 |
|
|
| 755 |
|
For a system with $f$ degrees of freedom the Liouville operator can be |
| 756 |
|
defined as, |
| 800 |
|
\label{introEq:Lp5} |
| 801 |
|
\end{align} |
| 802 |
|
Because $U_1(t)$ and $U_2(t)$ are unitary, $G(\Delta t)$ is also |
| 803 |
< |
unitary. Meaning an integrator based on this factorization will be |
| 803 |
> |
unitary. This means that an integrator based on this factorization will be |
| 804 |
|
reversible in time. |
| 805 |
|
|
| 806 |
|
As an example, consider the following decomposition of $L$: |
| 863 |
|
|
| 864 |
|
Based on the factorization from the previous section, |
| 865 |
|
Dullweber\emph{et al}.\cite{Dullweber1997}~ proposed a scheme for the |
| 866 |
< |
symplectic propagation of the rotation matrix, $\mathbf{A}$, as an |
| 866 |
> |
symplectic propagation of the rotation matrix, $\mathsf{A}$, as an |
| 867 |
|
alternative method for the integration of orientational degrees of |
| 868 |
|
freedom. The method starts with a straightforward splitting of the |
| 869 |
|
Liouville operator: |
| 870 |
|
\begin{align} |
| 871 |
|
iL_{\text{pos}} &= \dot{q}\frac{\partial}{\partial q} + |
| 872 |
< |
\mathbf{\dot{A}}\frac{\partial}{\partial \mathbf{A}} |
| 872 |
> |
\mathsf{\dot{A}}\frac{\partial}{\partial \mathsf{A}} |
| 873 |
|
\label{introEq:SR1a} \\% |
| 874 |
|
% |
| 875 |
|
iL_F &= F(q)\frac{\partial}{\partial p} |
| 876 |
|
\label{introEq:SR1b} \\% |
| 877 |
< |
iL_{\tau} &= \tau(\mathbf{A})\frac{\partial}{\partial \pi} |
| 877 |
> |
iL_{\tau} &= \tau(\mathsf{A})\frac{\partial}{\partial j} |
| 878 |
|
\label{introEq:SR1b} \\% |
| 879 |
|
\end{align} |
| 880 |
< |
Where $\tau(\mathbf{A})$ is the torque of the system |
| 881 |
< |
due to the configuration, and $\pi$ is the conjugate |
| 880 |
> |
Where $\tau(\mathsf{A})$ is the torque of the system |
| 881 |
> |
due to the configuration, and $j$ is the conjugate |
| 882 |
|
angular momenta of the system. The propagator, $G(\Delta t)$, becomes |
| 883 |
|
\begin{equation} |
| 884 |
|
G(\Delta t) = e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}} \, |
| 885 |
< |
e^{\frac{\Delta t}{2} \tau(\mathbf{A})\frac{\partial}{\partial \pi}} \, |
| 885 |
> |
e^{\frac{\Delta t}{2} \tau(\mathsf{A})\frac{\partial}{\partial j}} \, |
| 886 |
|
e^{\Delta t\,iL_{\text{pos}}} \, |
| 887 |
< |
e^{\frac{\Delta t}{2} \tau(\mathbf{A})\frac{\partial}{\partial \pi}} \, |
| 887 |
> |
e^{\frac{\Delta t}{2} \tau(\mathsf{A})\frac{\partial}{\partial j}} \, |
| 888 |
|
e^{\frac{\Delta t}{2} F(q)\frac{\partial}{\partial p}} |
| 889 |
|
\label{introEq:SR2} |
| 890 |
|
\end{equation} |
| 902 |
|
\mathcal{U}_x \biggl(\frac{\Delta t}{2}\biggr)\, |
| 903 |
|
\label{introEq:SR3} |
| 904 |
|
\end{equation} |
| 905 |
< |
Where $\mathcal{U}_j$ is a unitary rotation of $\mathbf{A}$ and |
| 906 |
< |
$\pi$ about each axis $j$. As all propagations are now |
| 905 |
> |
Where $\mathcal{U}_{\alpha}$ is a unitary rotation of $\mathsf{A}$ and |
| 906 |
> |
$j$ about each axis $\alpha$. As all propagations are now |
| 907 |
|
unitary and symplectic, the entire integration scheme is also |
| 908 |
|
symplectic and time reversible. |
| 909 |
|
|
| 910 |
|
\section{\label{introSec:layout}Dissertation Layout} |
| 911 |
|
|
| 912 |
< |
This dissertation is divided as follows:Ch.~\ref{chapt:RSA} |
| 912 |
> |
This dissertation is divided as follows: Ch.~\ref{chapt:RSA} |
| 913 |
|
presents the random sequential adsorption simulations of related |
| 914 |
|
pthalocyanines on a gold (111) surface. Ch.~\ref{chapt:oopse} |
| 915 |
< |
is about the writing of the molecular dynamics simulation package |
| 915 |
> |
is about our molecular dynamics simulation package |
| 916 |
|
{\sc oopse}. Ch.~\ref{chapt:lipid} regards the simulations of |
| 917 |
|
phospholipid bilayers using a mesoscale model. And lastly, |
| 918 |
|
Ch.~\ref{chapt:conclusion} concludes this dissertation with a |
| 924 |
|
study in applying Statistical Mechanics simulation techniques in order |
| 925 |
|
to obtain a simple model capable of explaining the results. My |
| 926 |
|
advisor, Dr. Gezelter, and I were approached by a colleague, |
| 927 |
< |
Dr. Lieberman, about possible explanations for the partial coverage of a |
| 928 |
< |
gold surface by a particular compound of hers. We suggested it might |
| 929 |
< |
be due to the statistical packing fraction of disks on a plane, and |
| 930 |
< |
set about to simulate this system. As the events in our model were |
| 931 |
< |
not dynamic in nature, a Monte Carlo method was employed. Here, if a |
| 932 |
< |
molecule landed on the surface without overlapping another, then its |
| 933 |
< |
landing was accepted. However, if there was overlap, the landing we |
| 934 |
< |
rejected and a new random landing location was chosen. This defined |
| 935 |
< |
our acceptance rules and allowed us to construct a Markov chain whose |
| 936 |
< |
limiting distribution was the surface coverage in which we were |
| 937 |
< |
interested. |
| 927 |
> |
Dr. Lieberman, about possible explanations for the partial coverage of |
| 928 |
> |
a gold surface by a particular compound synthesized in her group. We |
| 929 |
> |
suggested it might be due to the statistical packing fraction of disks |
| 930 |
> |
on a plane, and set about simulating this system. As the events in |
| 931 |
> |
our model were not dynamic in nature, a Monte Carlo method was |
| 932 |
> |
employed. Here, if a molecule landed on the surface without |
| 933 |
> |
overlapping with another, then its landing was accepted. However, if there |
| 934 |
> |
was overlap, the landing we rejected and a new random landing location |
| 935 |
> |
was chosen. This defined our acceptance rules and allowed us to |
| 936 |
> |
construct a Markov chain whose limiting distribution was the surface |
| 937 |
> |
coverage in which we were interested. |
| 938 |
|
|
| 939 |
|
The following chapter, about the simulation package {\sc oopse}, |
| 940 |
|
describes in detail the large body of scientific code that had to be |
| 941 |
|
written in order to study phospholipid bilayers. Although there are |
| 942 |
|
pre-existing molecular dynamic simulation packages available, none |
| 943 |
< |
were capable of implementing the models we were developing.{\sc oopse} |
| 944 |
< |
is a unique package capable of not only integrating the equations of |
| 945 |
< |
motion in Cartesian space, but is also able to integrate the |
| 946 |
< |
rotational motion of rigid bodies and dipoles. Add to this the |
| 947 |
< |
ability to perform calculations across parallel processors and a |
| 948 |
< |
flexible script syntax for creating systems, and {\sc oopse} becomes a |
| 949 |
< |
very powerful scientific instrument for the exploration of our model. |
| 943 |
> |
were capable of implementing the models we were developing. {\sc |
| 944 |
> |
oopse} is a unique package capable of not only integrating the |
| 945 |
> |
equations of motion in Cartesian space, but is also able to integrate |
| 946 |
> |
the rotational motion of rigid bodies and dipoles. It also has the |
| 947 |
> |
ability to perform calculations across distributed parallel processors |
| 948 |
> |
and contains a flexible script syntax for creating systems. {\sc |
| 949 |
> |
oopse} has become a very powerful scientific instrument for the |
| 950 |
> |
exploration of our model. |
| 951 |
|
|
| 952 |
< |
Bringing us to Ch.~\ref{chapt:lipid}. Using {\sc oopse}, I have been |
| 952 |
> |
In Ch.~\ref{chapt:lipid}, utilizing {\sc oopse}, I have been |
| 953 |
|
able to parameterize a mesoscale model for phospholipid simulations. |
| 954 |
|
This model retains information about solvent ordering around the |
| 955 |
|
bilayer, as well as information regarding the interaction of the |
| 959 |
|
provide the foundation for future exploration of bilayer phase |
| 960 |
|
behavior with this model. |
| 961 |
|
|
| 962 |
< |
Which leads into the last chapter, where I discuss future directions |
| 963 |
< |
for both{\sc oopse} and this mesoscale model. Additionally, I will |
| 964 |
< |
give a summary of results for this dissertation. |
| 962 |
> |
In the last chapter, I discuss future directions |
| 963 |
> |
for both {\sc oopse} and this mesoscale model. Additionally, I will |
| 964 |
> |
give a summary of the results found in this dissertation. |
| 965 |
|
|
| 966 |
|
|
