| 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}). |
| 339 |
> |
With the use of an $fix$, however, comes a discontinuity in the |
| 340 |
> |
potential energy curve (Fig.~\ref{fix}). To fix this discontinuity, |
| 341 |
> |
one calculates the potential energy at the $r_{\text{cut}}$, and add |
| 342 |
> |
that value to the potential. This causes the function to go smoothly |
| 343 |
> |
to zero at the cutoff radius. This ensures conservation of energy |
| 344 |
> |
when integrating the Newtonian equations of motion. |
| 345 |
> |
|
| 346 |
> |
The second main simplification used in this research is the Verlet |
| 347 |
> |
neighbor list. \cite{allen87:csl} In the Verlet method, one generates |
| 348 |
> |
a list of all neighbor atoms, $j$, surrounding atom $i$ within some |
| 349 |
> |
cutoff $r_{\text{list}}$, where $r_{\text{list}}>r_{\text{cut}}$. |
| 350 |
> |
This list is created the first time forces are evaluated, then on |
| 351 |
> |
subsequent force evaluations, pair calculations are only calculated |
| 352 |
> |
from the neighbor lists. The lists are updated if any given particle |
| 353 |
> |
in the system moves farther than $r_{\text{list}}-r_{\text{cut}}$, |
| 354 |
> |
giving rise to the possibility that a particle has left or joined a |
| 355 |
> |
neighbor list. |
| 356 |
|
|
| 357 |
+ |
\subsection{\label{introSec:MDintegrate} Integration of the equations of motion} |
| 358 |
|
|
| 359 |
+ |
A starting point for the discussion of molecular dynamics integrators |
| 360 |
+ |
is the Verlet algorithm. \cite{Frenkel1996} It begins with a Taylor |
| 361 |
+ |
expansion of position in time: |
| 362 |
+ |
\begin{equation} |
| 363 |
+ |
eq here |
| 364 |
+ |
\label{introEq:verletForward} |
| 365 |
+ |
\end{equation} |
| 366 |
+ |
As well as, |
| 367 |
+ |
\begin{equation} |
| 368 |
+ |
eq here |
| 369 |
+ |
\label{introEq:verletBack} |
| 370 |
+ |
\end{equation} |
| 371 |
+ |
Adding together Eq.~\ref{introEq:verletForward} and |
| 372 |
+ |
Eq.~\ref{introEq:verletBack} results in, |
| 373 |
+ |
\begin{equation} |
| 374 |
+ |
eq here |
| 375 |
+ |
\label{introEq:verletSum} |
| 376 |
+ |
\end{equation} |
| 377 |
+ |
Or equivalently, |
| 378 |
+ |
\begin{equation} |
| 379 |
+ |
eq here |
| 380 |
+ |
\label{introEq:verletFinal} |
| 381 |
+ |
\end{equation} |
| 382 |
+ |
Which contains an error in the estimate of the new positions on the |
| 383 |
+ |
order of $\Delta t^4$. |
| 384 |
+ |
|
| 385 |
+ |
In practice, however, the simulations in this research were integrated |
| 386 |
+ |
with a velocity reformulation of teh Verlet method. \cite{allen87:csl} |
| 387 |
+ |
\begin{equation} |
| 388 |
+ |
eq here |
| 389 |
+ |
\label{introEq:MDvelVerletPos} |
| 390 |
+ |
\end{equation} |
| 391 |
+ |
\begin{equation} |
| 392 |
+ |
eq here |
| 393 |
+ |
\label{introEq:MDvelVerletVel} |
| 394 |
+ |
\end{equation} |
| 395 |
+ |
The original Verlet algorithm can be regained by substituting the |
| 396 |
+ |
velocity back into Eq.~\ref{introEq:MDvelVerletPos}. The Verlet |
| 397 |
+ |
formulations are chosen in this research because the algorithms have |
| 398 |
+ |
very little long term drift in energy conservation. Energy |
| 399 |
+ |
conservation in a molecular dynamics simulation is of extreme |
| 400 |
+ |
importance, as it is a measure of how closely one is following the |
| 401 |
+ |
``true'' trajectory wtih the finite integration scheme. An exact |
| 402 |
+ |
solution to the integration will conserve area in phase space, as well |
| 403 |
+ |
as be reversible in time, that is, the trajectory integrated forward |
| 404 |
+ |
or backwards will exactly match itself. Having a finite algorithm |
| 405 |
+ |
that both conserves area in phase space and is time reversible, |
| 406 |
+ |
therefore increases, but does not guarantee the ``correctness'' or the |
| 407 |
+ |
integrated trajectory. |
| 408 |
+ |
|
| 409 |
+ |
It can be shown, \cite{Frenkel1996} that although the Verlet algorithm |
| 410 |
+ |
does not rigorously preserve the actual Hamiltonian, it does preserve |
| 411 |
+ |
a pseudo-Hamiltonian which shadows the real one in phase space. This |
| 412 |
+ |
pseudo-Hamiltonian is proveably area-conserving as well as time |
| 413 |
+ |
reversible. The fact that it shadows the true Hamiltonian in phase |
| 414 |
+ |
space is acceptable in actual simulations as one is interested in the |
| 415 |
+ |
ensemble average of the observable being measured. From the ergodic |
| 416 |
+ |
hypothesis (Sec.~\ref{introSec:StatThermo}), it is known that the time |
| 417 |
+ |
average will match the ensemble average, therefore two similar |
| 418 |
+ |
trajectories in phase space should give matching statistical averages. |
| 419 |
+ |
|
| 420 |
+ |
\subsection{\label{introSec:MDfurther}Further Considerations} |
| 421 |
+ |
In the simulations presented in this research, a few additional |
| 422 |
+ |
parameters are needed to describe the motions. The simulations |
| 423 |
+ |
involving water and phospholipids in Chapt.~\ref{chaptLipids} are |
| 424 |
+ |
required to integrate the equations of motions for dipoles on atoms. |
| 425 |
+ |
This involves an additional three parameters be specified for each |
| 426 |
+ |
dipole atom: $\phi$, $\theta$, and $\psi$. These three angles are |
| 427 |
+ |
taken to be the Euler angles, where $\phi$ is a rotation about the |
| 428 |
+ |
$z$-axis, and $\theta$ is a rotation about the new $x$-axis, and |
| 429 |
+ |
$\psi$ is a final rotation about the new $z$-axis (see |
| 430 |
+ |
Fig.~\ref{introFig:euleerAngles}). This sequence of rotations can be |
| 431 |
+ |
accumulated into a single $3 \times 3$ matrix $\mathbf{A}$ |
| 432 |
+ |
defined as follows: |
| 433 |
+ |
\begin{equation} |
| 434 |
+ |
eq here |
| 435 |
+ |
\label{introEq:EulerRotMat} |
| 436 |
+ |
\end{equation} |
| 437 |
+ |
|
| 438 |
+ |
The equations of motion for Euler angles can be written down as |
| 439 |
+ |
\cite{allen87:csl} |
| 440 |
+ |
\begin{equation} |
| 441 |
+ |
eq here |
| 442 |
+ |
\label{introEq:MDeuleeerPsi} |
| 443 |
+ |
\end{equation} |
| 444 |
+ |
Where $\omega^s_i$ is the angular velocity in the lab space frame |
| 445 |
+ |
along cartesian coordinate $i$. However, a difficulty arises when |
| 446 |
+ |
attempting to integrate Eq.~\ref{introEq:MDeulerPhi} and |
| 447 |
+ |
Eq.~\ref{introEq:MDeulerPsi}. The $\frac{1}{\sin \theta}$ present in |
| 448 |
+ |
both equations means there is a non-physical instability present when |
| 449 |
+ |
$\theta$ is 0 or $\pi$. |
| 450 |
+ |
|
| 451 |
+ |
To correct for this, the simulations integrate the rotation matrix, |
| 452 |
+ |
$\mathbf{A}$, directly, thus avoiding the instability. |
| 453 |
+ |
This method was proposed by Dullwebber |
| 454 |
+ |
\emph{et. al.}\cite{Dullwebber:1997}, and is presented in |
| 455 |
+ |
Sec.~\ref{introSec:MDsymplecticRot}. |
| 456 |
+ |
|
| 457 |
+ |
\subsubsection{\label{introSec:MDliouville}Liouville Propagator} |
| 458 |
+ |
|
| 459 |
+ |
|
| 460 |
|
\section{\label{introSec:chapterLayout}Chapter Layout} |
| 461 |
|
|
| 462 |
|
\subsection{\label{introSec:RSA}Random Sequential Adsorption} |
| 463 |
|
|
| 464 |
|
\subsection{\label{introSec:OOPSE}The OOPSE Simulation Package} |
| 465 |
|
|
| 466 |
< |
\subsection{\label{introSec:bilayers}A Mesoscale Model for Phospholipid Bilayers} |
| 466 |
> |
\subsection{\label{introSec:bilayers}A Mesoscale Model for |
| 467 |
> |
Phospholipid Bilayers} |
