trunk/chrisDissertation/Introduction.tex

\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}

In order to advance our understanding of natural chemical and physical
processes, researchers develop explanations for events observed
experimentally. These hypotheses, supported by a body of corroborating
observations, can develop into accepted theories for how these
processes occur. This validation process, as well as testing the
limits of these theories, is essential in developing a firm foundation
for our knowledge. Theories involving molecular scale systems cause
particular difficulties in this process because their size and often
fast motions make them difficult to observe experimentally. One useful
tool for addressing these difficulties is computer simulation.
 
In computer simulations, we can develop models from either the theory
or experimental knowledge and then test them in a controlled
environment. Work done in this manner allows us to further refine
theories, more accurately represent what is happening in experimental
observations, and even make predictions about what one will see in
experiments. Thus, computer simulations of molecular systems act as a
bridge between theory and experiment.

Depending on the system of interest, there are a variety of different
computational techniques that can used to test and gather information
from the developed models. In the study of classical systems, the two
most commonly used techniques are Monte Carlo and molecular
dynamics. Both of these methods operate by calculating interactions
between particles of our model systems; however, the progression of
the simulation under the different techniques is vastly
different. Monte Carlo operates through random configuration changes
that that follow rules adhering to a specific statistical mechanics
ensemble, while molecular dynamics is chiefly concerned with solving
the classic equations of motion to move between configurations within
an ensemble. Thermodynamic properties can be calculated from both
techniques; but because of the random nature of Monte Carlo, only
molecular dynamics can be used to investigate dynamical
quantities. The research presented in the following chapters utilized
molecular dynamics near exclusively, so we will present a general
introduction to molecular dynamics and not Monte Carlo. There are
several resources available for those desiring a more in-depth
presentation of either of these
techniques.\cite{Allen87,Frenkel02,Leach01}

\section{\label{sec:MolecularDynamics}Molecular Dynamics}

As stated above, in molecular dynamics we focus on evolving
configurations of molecules over time. In order to use this as a tool
for understanding experiments and testing theories, we want the
configuration to evolve in a manner that mimics real molecular
systems. To do this, we start with clarifying what we know about a
given configuration of particles at time $t_1$, basically that each
particle in the configuration has a position ($q$) and velocity
($\dot{q}$). We now want to know what the configuration will be at
time $t_2$. To find out, we need the classical equations of
motion, and one useful formulation of them is the Lagrangian form.

The Lagrangian ($L$) is a function of the positions and velocities that
takes the form,
\begin{equation}
L = K - V,
\label{eq:lagrangian}
\end{equation}
where $K$ is the kinetic energy and $V$ is the potential energy. We
can use Hamilton's principle, which states that the integral of the
Lagrangian over time has a stationary value for the correct path of
motion, to say that the variation of the integral of the Lagrangian
over time is zero,\cite{Tolman38}
\begin{equation}
\delta\int_{t_1}^{t_2}L(q,\dot{q})dt = 0.
\end{equation}
This can be written as a summation of integrals to give
\begin{equation}
\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left(\frac{\partial L}{\partial q_i}\delta q_i
        + \frac{\partial L}{\partial \dot{q}_i}\delta \dot{q}_i\right)dt = 0.
\end{equation}
Using fact that $\dot q$ is the derivative with respect to time of $q$
and integrating the second partial derivative in the parenthesis by
parts, this equation simplifies to
\begin{equation}
\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left(
        \frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i}
        - \frac{\partial L}{\partial q_i}\right)\delta {q}_i dt = 0,
\end{equation}
and the above equation is only valid for
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial \dot{q}_i}
        - \frac{\partial L}{\partial q_i} = 0\quad\quad(i=1,2,\dots,3N).
\label{eq:formulation}
\end{equation}
To obtain the classical equations of motion for the particles, we can
substitute equation (\ref{eq:lagrangian}) into the above equation with
$m\dot{\mathbf{r}}^2/2$ for the kinetic energy, giving
\begin{equation}
\frac{d}{dt}(m\dot{\mathbf{r}})+\frac{dV}{d\mathbf{r}}=0,
\end{equation}
or more recognizably,
\begin{equation}
\mathbf{f} = m\mathbf{a},
\end{equation}
where $\mathbf{f} = -dV/d\mathbf{r}$ and $\mathbf{a} =
d^2\mathbf{r}/dt^2$. This Lagrangian formulation shown in equation
(\ref{eq:formulation}) is generalized, and it can be used to determine
equations of motion in forms outside of the typical Cartesian case
shown here.

\subsection{\label{sec:Verlet}Verlet Integration}

In order to perform molecular dynamics, we need an algorithm that
integrates the equations of motion described above. Ideal algorithms
are both simple in implementation and highly accurate. There have been
a large number of algorithms developed for this purpose; however, for
reasons discussed below, we are going to focus on the Verlet class of
integrators.\cite{Gear66,Beeman76,Berendsen86,Allen87,Verlet67,Swope82}

In Verlet's original study of computer ``experiments'', he directly
integrated the Newtonian second order differential equation of motion,
\begin{equation}
m\frac{d^2\mathbf{r}_i}{dt^2} = \sum_{j\ne i}\mathbf{f}(r_{ij}),
\end{equation}
with the following simple algorithm:
\begin{equation}
\mathbf{r}_i(t+\delta t) = -\mathbf{r}_i(t-\delta t) + 2\mathbf{r}_i(t)
        + \sum_{j\ne i}\mathbf{f}(r_{ij}(t))\delta t^2,
\label{eq:verlet}
\end{equation}
where $\delta t$ is the time step of integration.\cite{Verlet67} It is
interesting to note that equation (\ref{eq:verlet}) does not include
velocities, and this makes sense since they are not present in the
differential equation. The velocities are necessary for the
calculation of the kinetic energy, and they can be calculated at each
time step with the following equation:
\begin{equation}
\mathbf{v}_i(t) = \frac{\mathbf{r}_i(t+\delta t)-\mathbf{r}_i(t-\delta t)}
                       {2\delta t}.
\end{equation}

Like the equation of motion it solves, the Verlet algorithm has the
beneficial property of being time-reversible, meaning that you can
integrate the configuration forward and then backward and end up at
the original configuration. Some other methods for integration, like
predictor-corrector methods, lack this property in that they require
higher order information that is discarded after integrating
steps. Another interesting property of this algorithm is that it is
symplectic, meaning that it preserves area in phase-space. Symplectic
algorithms system stays in the region of phase-space dictated by the
ensemble and enjoy greater long-time energy
conservations.\cite{Frenkel02}

While the error in the positions calculated using the Verlet algorithm
is small ($\mathcal{O}(\delta t^4)$), the error in the velocities is
quite a bit larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope
{\it et al.}  developed a corrected for of this algorithm, the
`velocity Verlet' algorithm, which improves the error of the velocity
calculation and thus the energy conservation.\cite{Swope82} This
algorithm involves a full step of the positions,
\begin{equation}
\mathbf{r}(t+\delta t) = \mathbf{r}(t) + \delta t\mathbf{v}(t)
                                + \frac{1}{2}\delta t^2\mathbf{a}(t),
\end{equation}
and a half step of the velocities,
\begin{equation}
\mathbf{v}\left(t+\frac{1}{2}\delta t\right) = \mathbf{v}(t) 
                                        + \frac{1}{2}\delta t\mathbf{a}(t).
\end{equation}
After calculating new forces, the velocities can be updated to a full
step,
\begin{equation}
\mathbf{v}(t+\delta t) = \mathbf{v}\left(t+\frac{1}{2}\delta t\right) 
                                + \frac{1}{2}\delta t\mathbf{a}(t+\delta t).
\end{equation}
By integrating in this manner, the error in the velocities reduces to
$\mathcal{O}(\delta t^3)$. It should be noted that the error in the
positions increases to $\mathcal{O}(\delta t^3)$, but the resulting
improvement in the energies coupled with the maintained simplicity,
time-reversibility, and symplectic nature make it an improvement over
the original form.

\subsection{\label{sec:IntroIntegrate}Rigid Body Motion}

Rigid bodies are non-spherical particles or collections of particles
(e.g. $\mbox{C}_{60}$) that have a constant internal potential and
move collectively.\cite{Goldstein01} Discounting iterative constraint
procedures like {\sc shake} and {\sc rattle} for approximating rigid
bodies, they are not included in most simulation packages because of
the algorithmic complexity involved in propagating orientational
degrees of freedom.\cite{Ryckaert77,Andersen83,Krautler01} Integrators
which propagate orientational motion with an acceptable level of
energy conservation for molecular dynamics are relatively new
inventions.

Moving a rigid body involves determination of both the force and
torque applied by the surroundings, which directly affect the
translational and rotational motion in turn. In order to accumulate
the total force on a rigid body, the external forces and torques must
first be calculated for all the internal particles. The total force on
the rigid body is simply the sum of these external forces.
Accumulation of the total torque on the rigid body is more complex
than the force because the torque is applied to the center of mass of
the rigid body. The space-fixed torque on rigid body $i$ is
\begin{equation}
\boldsymbol{\tau}_i=
        \sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} 
        + \boldsymbol{\tau}_{ia}\biggr],
\label{eq:torqueAccumulate}
\end{equation}
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
position of the center of mass respectively, while $\mathbf{f}_{ia}$,
$\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on,
position of, and torque on the component particles.

The summation of the total torque is done in the body fixed axis. In
order to move between the space fixed and body fixed coordinate axes,
parameters describing the orientation must be maintained for each
rigid body. At a minimum, the rotation matrix ($\mathsf{A}$) can be
described by the three Euler angles ($\phi, \theta,$ and $\psi$),
where the elements of $\mathsf{A}$ are composed of trigonometric
operations involving $\phi, \theta,$ and $\psi$.\cite{Goldstein01}
Direct propagation of the Euler angles has a known $1/\sin\theta$
divergence in the equations of motion for $\phi$ and $\psi$, leading
to numerical instabilities any time one of the directional atoms or
rigid bodies has an orientation near $\theta=0$ or
$\theta=\pi$.\cite{Allen87} One of the most practical ways to avoid
this ``gimbal point'' is to switch to another angular set defining the
orientation of the rigid body near this point.\cite{Barojas73} This
procedure results in additional book-keeping and increased algorithm
complexity. In the search for more elegant alternative methods, Evans
proposed the use of quaternions to describe and propagate
orientational motion.\cite{Evans77}

The quaternion method for integration involves a four dimensional
representation of the orientation of a rigid
body.\cite{Evans77,Evans77b,Allen87} Thus, the elements of
$\mathsf{A}$ can be expressed as arithmetic operations involving the
four quaternions ($q_w, q_x, q_y,$ and $q_z$),
\begin{equation}
\mathsf{A} = \left( \begin{array}{l@{\quad}l@{\quad}l}
q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \\
2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \\
2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2 \\
\end{array}\right).
\end{equation}
Integration of the equations of motion involves a series of arithmetic
operations involving the quaternions and angular momenta and leads to
performance enhancements over Euler angles, particularly for very
small systems.\cite{Evans77} This integration methods works well for
propagating orientational motion in the canonical ensemble ($NVT$);
however, energy conservation concerns arise when using the simple
quaternion technique under the microcanonical ($NVE$) ensemble.  An
earlier implementation of {\sc oopse} utilized quaternions for
propagation of rotational motion; however, a detailed investigation
showed that they resulted in a steady drift in the total energy,
something that has been observed by Kol {\it et al.} (also see
section~\ref{sec:waterSimMethods}).\cite{Kol97}

Because of the outlined issues involving integration of the
orientational motion using both Euler angles and quaternions, we
decided to focus on a relatively new scheme that propagates the entire
nine parameter rotation matrix. This techniques is a velocity-Verlet
version of the symplectic splitting method proposed by Dullweber,
Leimkuhler and McLachlan ({\sc dlm}).\cite{Dullweber97} When there are
no directional atoms or rigid bodies present in the simulation, this
integrator becomes the standard velocity-Verlet integrator which is
known to effectively sample the microcanonical ($NVE$)
ensemble.\cite{Frenkel02}

The key aspect of the integration method proposed by Dullweber
\emph{et al.} is that the entire $3 \times 3$ rotation matrix is
propagated from one time step to the next. In the past, this would not
have been as feasible, since the rotation matrix for a single body has
nine elements compared with the more memory-efficient methods (using
three Euler angles or four quaternions).  Computer memory has become
much less costly in recent years, and this can be translated into
substantial benefits in energy conservation.

The integration of the equations of motion is carried out in a
velocity-Verlet style two-part algorithm.\cite{Swope82} The first-part
({\tt moveA}) consists of a half-step ($t + \delta t/2$) of the linear
velocity (${\bf v}$) and angular momenta ({\bf j}) and a full-step ($t
+ \delta t$) of the positions ({\bf r}) and rotation matrix,
\begin{equation*}
{\tt moveA} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l}
{\bf v}\left(t + \delta t / 2\right) & {\bf v}(t) 
        + \left( {\bf f}(t) / m \right)(\delta t/2), \\
%
{\bf r}(t + \delta t) & {\bf r}(t) 
        + {\bf v}\left(t + \delta t / 2 \right)\delta t, \\
%
{\bf j}\left(t + \delta t / 2 \right) & {\bf j}(t) 
        + \boldsymbol{\tau}^b(t)(\delta t/2), \\
%
\mathsf{A}(t + \delta t) & \mathrm{rotate}\left( {\bf j}
        (t + \delta t / 2)\delta t \cdot 
                \overleftrightarrow{\mathsf{I}}^{-1} \right),
\end{array}\right.
\end{equation*}
where $\overleftrightarrow{\mathsf{I}}^{-1}$ is the inverse of the
moment of inertia tensor. The $\mathrm{rotate}$ function is the
product of rotations about the three body-fixed axes,
\begin{equation}
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y /
2) \cdot \mathsf{G}_x(a_x /2),
\label{eq:dlmTrot}
\end{equation}
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates
both the rotation matrix ($\mathsf{A}$) and the body-fixed angular
momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis
$\alpha$,
\begin{equation}
\mathsf{G}_\alpha( \theta ) = \left\{
\begin{array}{l@{\quad\leftarrow\quad}l}
\mathsf{A}(t) & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^\textrm{T},\\
{\bf j}(t) & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0).
\end{array}
\right.
\end{equation}
$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
rotation matrix.  For example, in the small-angle limit, the rotation
matrix around the body-fixed x-axis can be approximated as
\begin{equation}
\mathsf{R}_x(\theta) \approx \left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+\theta^2 / 4} \\
0 & \frac{\theta}{1+\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}
\end{array}
\right).
\end{equation}
The remaining rotations follow in a straightforward manner. As seen
from the form of equation~(\ref{eq:dlmTrot}), the {\sc dlm} method
uses a Trotter factorization of the orientational
propagator.\cite{Trotter59} This has three effects:
\begin{enumerate}
\item the integrator is area-preserving in phase space (i.e. it is
{\it symplectic}),
\item the integrator is time-{\it reversible}, and
\item the error for a single time step is of order 
$\mathcal{O}\left(\delta t^4\right)$ for time steps of length $\delta t$.
\end{enumerate}

After the initial half-step ({\tt moveA}), the forces and torques are
evaluated for all of the particles. Once completed, the velocities can
be advanced to complete the second-half of the two-part algorithm
({\tt moveB}), resulting an a completed full step of both the
positions and momenta,
\begin{equation*}
{\tt moveB} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l}
{\bf v}\left(t + \delta t \right) &
        {\bf v}\left(t + \delta t / 2 \right) 
        + \left({\bf f}(t + \delta t) / m \right)(\delta t/2), \\
%
{\bf j}\left(t + \delta t \right) &
        {\bf j}\left(t + \delta t / 2 \right) 
        + \boldsymbol{\tau}^b(t + \delta t)(\delta t/2).
\end{array}\right.
\end{equation*}

The matrix rotations used in the {\sc dlm} method end up being more
costly computationally than the simpler arithmetic quaternion
propagation. With the same time step, a 1024-molecule water simulation
incurs approximately a 10\% increase in computation time using the
{\sc dlm} method in place of quaternions. This cost is more than
justified when comparing the energy conservation achieved by the two
methods. Figure \ref{fig:quatdlm} provides a comparative analysis of
the {\sc dlm} method versus the traditional quaternion scheme.

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/dlmVsQuat.pdf}
\caption[Energy conservation analysis of the {\sc dlm} and quaternion 
integration methods]{Analysis of the energy conservation of the {\sc
dlm} and quaternion integration methods.  $\delta \mathrm{E}_1$ is the
linear drift in energy over time and $\delta \mathrm{E}_0$ is the
standard deviation of energy fluctuations around this drift.  All
simulations were of a 1024 SSD water system at 298 K starting from the
same initial configuration. Note that the {\sc dlm} method provides
more than an order-of-magnitude improvement in both the energy drift
and the size of the energy fluctuations when compared with the
quaternion method at any given time step.  At time steps larger than 4
fs, the quaternion scheme resulted in rapidly rising energies which
eventually lead to simulation failure.  Using the {\sc dlm} method,
time steps up to 8 fs can be taken before this behavior is evident.}
\label{fig:quatdlm}
\end{figure}

In figure \ref{fig:quatdlm}, $\delta \mbox{E}_1$ is a measure of the
linear energy drift in units of $\mbox{kcal mol}^{-1}$ per particle
over a nanosecond of simulation time, and $\delta \mbox{E}_0$ is the
standard deviation of the energy fluctuations in units of $\mbox{kcal
mol}^{-1}$ per particle. In the top plot, it is apparent that the
energy drift is reduced by a significant amount (2 to 3 orders of
magnitude improvement at all tested time steps) by choosing the {\sc
dlm} method over the simple non-symplectic quaternion integration
method.  In addition to this improvement in energy drift, the
fluctuations in the total energy are also dampened by 1 to 2 orders of
magnitude by utilizing the {\sc dlm} method.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/compCost.pdf}
\caption[Energy drift as a function of required simulation run 
time]{Energy drift as a function of required simulation run time.
$\delta \mathrm{E}_1$ is the linear drift in energy over time.
Simulations were performed on a single 2.5 GHz Pentium 4
processor. Simulation time comparisons can be made by tracing
horizontally from one curve to the other. For example, a simulation
that takes 24 hours using the {\sc dlm} method will take roughly
210 hours using the simple quaternion method if the same degree of
energy conservation is desired.}
\label{fig:cpuCost}
\end{figure}
Although the {\sc dlm} method is more computationally expensive than
the traditional quaternion scheme for propagating of a time step,
consideration of the computational cost for a long simulation with a
particular level of energy conservation is in order.  A plot of energy
drift versus computational cost was generated
(Fig.~\ref{fig:cpuCost}). This figure provides an estimate of the CPU
time required under the two integration schemes for 1 nanosecond of
simulation time for the model 1024-molecule system.  By choosing a
desired energy drift value it is possible to determine the CPU time
required for both methods. If a $\delta \mbox{E}_1$ of
0.001~kcal~mol$^{-1}$ per particle is desired, a nanosecond of
simulation time will require ~19 hours of CPU time with the {\sc dlm}
integrator, while the quaternion scheme will require ~154 hours of CPU
time. This demonstrates the computational advantage of the integration
scheme utilized in {\sc oopse}.

\section{Accumulating Interactions}

In the force calculation between {\tt moveA} and {\tt moveB} mentioned
in section \ref{sec:IntroIntegrate}, we need to accumulate the
potential and forces (and torques if the particle is a rigid body or
multipole) on each particle from their surroundings. This can quickly
become a cumbersome task for large systems since the number of pair
interactions increases by $\mathcal{O}(N(N-1)/2)$ if you accumulate
interactions between all particles in the system, note the utilization
of Newton's third law to reduce the interaction number from
$\mathcal{O}(N^2)$. The case of periodic boundary conditions further
complicates matters by turning the finite system into an infinitely
repeating one. Fortunately, we can reduce the scale of this problem by
using spherical cutoff methods.

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/sphericalCut.pdf}
\caption{When using a spherical cutoff, only particles within a chosen
cutoff radius distance, $R_\textrm{c}$, of the central particle are
included in the pairwise summation. This reduces a problem that scales
by $\sim\mathcal{O}(N^2)$ to one that scales by $\sim\mathcal{O}(N)$.}
\label{fig:sphereCut}
\end{figure}
With spherical cutoffs, rather than accumulating the full set of
interactions between all particles in the simulation, we only
explicitly consider interactions between local particles out to a
specified cutoff radius distance, $R_\textrm{c}$, (see figure
\ref{fig:sphereCut}). This reduces the scaling of the interaction to
$\mathcal{O}(N\cdot\textrm{c})$, where `c' is a value that depends on
the size of $R_\textrm{c}$ (c $\approx R_\textrm{c}^3$). Determination
of which particles are within the cutoff is also an issue, because
this process requires a full loop over all $N(N-1)/2$ pairs. To reduce
the this expense, we can use neighbor lists.\cite{Verlet67,Thompson83}
With neighbor lists, we have a second list of particles built from a
list radius $R_\textrm{l}$, which is larger than $R_\textrm{c}$. Once
any of the particles within $R_\textrm{l}$ move a distance of
$R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we rebuild the
list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an appropriate
skin thickness, these updates are only performed every $\sim$20
time steps, significantly reducing the time spent on pair-list
bookkeeping operations.\cite{Allen87} If these neighbor lists are
utilized, it is important that these list updates occur
regularly. Incorrect application of this technique leads to
non-physical dynamics, such as the ``flying block of ice'' behavior
for which improper neighbor list handling was identified a one of the
possible causes.\cite{Harvey98,Sagui99}

\subsection{Correcting Cutoff Discontinuities}
\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/ljCutoffSquare.pdf}
\caption{The common methods to smooth the potential discontinuity 
introduced when using a cutoff include a shifted potential, a shifted
force, and a switching function. The shifted potential and shifted
force both lift the whole potential so that it zeroes at
$R_\textrm{c}$, thereby reducing the strength of the interaction. The
(cubic) switching function only alters the potential in the switching
region in order to smooth out the discontinuity.}
\label{fig:ljCutoff}
\end{figure}
As particle move in and out of $R_\textrm{c}$, there will be sudden
discontinuous jumps in the potential (and forces) due to their
appearance and disappearance. In order to prevent heating and poor
energy conservation in the simulations, this discontinuity needs to be
smoothed out. There are several ways to modify the function so that it
crosses $R_\textrm{c}$ in a continuous fashion, and the easiest
methods is shifting the potential. To shift the potential, we simply
subtract out the value we calculate at $R_\textrm{c}$ from the whole
potential. For the shifted form of the Lennard-Jones potential is
\begin{equation}
V_\textrm{SLJ} = \left\{\begin{array}{l@{\quad\quad}l}
        V_\textrm{LJ}(r_{ij}) - V_\textrm{LJ}(R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\
        0 & r_{ij} \geqslant R_\textrm{c},
\end{array}\right.
\end{equation}
where
\begin{equation}
V_\textrm{LJ} = 4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} -
                                \left(\frac{\sigma}{r_{ij}}\right)^6\right].
\end{equation}
In figure \ref{fig:ljCutoff}, the shifted form of the potential
reaches zero at the cutoff radius at the expense of the correct
magnitude below $R_\textrm{c}$. This correction method also does
nothing to correct the discontinuity in the forces. We can account for
this force discontinuity by shifting in the by applying the shifting
in the forces rather than just the potential via
\begin{equation}
V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l}
        V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) - 
                \left(\frac{d V(r_{ij})}{d r_{ij}}\right)_{r_{ij}=R_\textrm{c}}
                (r_{ij} - R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\
        0 & r_{ij} \geqslant R_\textrm{c}.
\end{array}\right.
\end{equation}
The forces are continuous with this potential; however, the inclusion
of the derivative term distorts the potential even further than the
simple shifting as shown in figure \ref{fig:ljCutoff}. The method for
correcting the discontinuity which results in the smallest
perturbation in both the potential and the forces is the use of a
switching function. The cubic switching function,
\begin{equation}
S(r) = \left\{\begin{array}{l@{\quad\quad}l}
        1 & r_{ij} \leqslant R_\textrm{sw}, \\
        \frac{(R_\textrm{c} + 2r_{ij} - 3R_\textrm{sw})
                (R_\textrm{c} - r_{ij})^2}{(R_\textrm{c} - R_\textrm{sw})^3} 
                & R_\textrm{sw} < r_{ij} \leqslant R_\textrm{c}, \\
        0 & r_{ij} > R_\textrm{c},
        \end{array}\right.
\end{equation}
is sufficient to smooth the potential (again, see figure
\ref{fig:ljCutoff}) and the forces by only perturbing the potential in
the switching region. If smooth second derivatives are required, a
higher order polynomial switching function (i.e. fifth order
polynomial) can be used.\cite{Andrea83,Leach01} It should be noted
that the higher the order of the polynomial, the more abrupt the
switching transition.

\subsection{\label{sec:LJCorrections}Long-Range Interaction Corrections}

While a good approximation, accumulating interaction only from the
nearby particles can lead to some issues, because the further away
surrounding particles do still have a small effect. For instance,
while the strength of the Lennard-Jones interaction is quite weak at
$R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding all of
the attractive interaction that extends out to extremely long
distances. Thus, the potential is a little too high and the pressure
on the central particle from the surroundings is a little too low. For
homogeneous Lennard-Jones systems, we can correct for this neglect by
assuming a uniform density and integrating the missing part,
\begin{equation}
V_\textrm{full}(r_{ij}) \approx V_\textrm{c} 
                + 2\pi N\rho\int^\infty_{R_\textrm{c}}r^2V_\textrm{LJ}(r)dr,
\end{equation} 
where $V_\textrm{c}$ is the truncated Lennard-Jones
potential.\cite{Allen87} Like the potential, other Lennard-Jones
properties can be corrected by integration over the relevant
function. Note that with heterogeneous systems, this correction begins
to break down because the density is no longer uniform.

Correcting long-range electrostatic interactions is a topic of great
importance in the field of molecular simulations. There have been
several techniques developed to address this issue, like the Ewald
summation and the reaction field technique. An in-depth analysis of
this problem, as well as useful corrective techniques, is presented in
chapter \ref{chap:electrostatics}.

\subsection{Periodic Boundary Conditions}

In typical molecular dynamics simulations there are no restrictions
placed on the motion of particles outside of what the inter-particle
interactions dictate. This means that if a particle collides with
other particles, it is free to move away from the site of the
collision. If we consider the entire system as a collection of
particles, they are not confined by walls of the ``simulation box''
and can freely move away from the other particles. With no boundary
considerations, particles moving outside of the simulation box
enter a vacuum. This is correct behavior for cluster simulations in a
vacuum; however, if we want to simulate bulk or spatially infinite
systems, we need to use periodic boundary conditions.

\begin{figure}
\centering
\includegraphics[width=4.5in]{./figures/periodicImage.pdf}
\caption{With periodic boundary conditions imposed, when particles 
move out of one side the simulation box, they wrap back in the
opposite side. In this manner, a finite system of particles behaves as
an infinite system.}
\label{fig:periodicImage}
\end{figure}
In periodic boundary conditions, as a particle moves outside one wall
of the simulation box, the coordinates are remapped such that the
particle enters the opposing side of the box. This process is easy to
visualize in two dimensions as shown in figure \ref{fig:periodicImage}
and can occur in three dimensions, though it is not as easy to
visualize. Remapping the actual coordinates of the particles can be
problematic in that the we are restricting the distance a particle can
move from it's point of origin to a diagonal of the simulation
box. Thus, even though we are not confining the system with hard
walls, we are confining the particle coordinates to a particular
region in space. To avoid this ``soft'' confinement, it is common
practice to allow the particle coordinates to expand in an
unrestricted fashion while calculating interactions using a wrapped
set of ``minimum image'' coordinates. These minimum image coordinates
need not be stored because they are easily calculated on the fly when
determining particle distances.

\section{Calculating Properties}

In order to use simulations to model experimental process and evaluate
theories, we need to be able to extract properties from the
results. In experiments, we can measure thermodynamic properties such
as the pressure and free energy. In computer simulations, we can
calculate properties from the motion and configuration of particles in
the system and make connections between these properties and the
experimental thermodynamic properties through statistical mechanics.

The work presented in the later chapters use the canonical ($NVT$),
isobaric-isothermal ($NPT$), and microcanonical ($NVE$) statistical
mechanics ensembles. The different ensembles lend themselves to more
effectively calculating specific properties. For instance, if we
concerned ourselves with the calculation of dynamic properties, which
are dependent upon the motion of the particles, it is better to choose
an ensemble that does not add system motions to keep properties like
the temperature or pressure constant. In this case, the $NVE$ ensemble
would be the most appropriate choice. In chapter \ref{chap:ice}, we
discuss calculating free energies, which are non-mechanical
thermodynamic properties, and these calculations also depend on the
chosen ensemble.\cite{Allen87} The Helmholtz free energy ($A$) depends
on the $NVT$ partition function ($Q_{NVT}$),
\begin{equation}
A = -k_\textrm{B}T\ln Q_{NVT},
\end{equation}
while the Gibbs free energy ($G$) depends on the $NPT$ partition
function ($Q_{NPT}$),
\begin{equation}
G = -k_\textrm{B}T\ln Q_{NPT}.  
\end{equation}
It is also useful to note that the conserved quantities of the $NVT$
and $NPT$ ensembles are related to the Helmholtz and Gibbs free
energies respectively.\cite{Melchionna93}

Integrating the equations of motion is a simple method to obtain a
sequence of configurations that sample the chosen ensemble. For each
of these configurations, we can calculate an instantaneous value for a
chosen property like the density in the $NPT$ ensemble, where the
volume is allowed to fluctuate. The density for the simulation is
calculated from an average over the densities for the individual
configurations. Being that the configurations from the integration
process are related to one another by the time evolution of the
interactions, this average is technically a time average. In
calculating thermodynamic properties, we would actually prefer an
ensemble average that is representative of all accessible states of
the system. We can calculate thermodynamic properties from the time
average by taking advantage of the ergodic hypothesis, which states
that over a long period of time, the time average and the ensemble
average are the same.

In addition to the average, the fluctuations of that particular
property can be determined via the standard deviation. Fluctuations
are useful for measuring various thermodynamic properties in computer
simulations. In section \ref{sec:t5peThermo}, we use fluctuations in
properties like the enthalpy and volume to calculate various
thermodynamic properties, such as the constant pressure heat capacity
and the isothermal compressibility.

\section{Application of the Techniques}

In the following chapters, the above techniques are all utilized in
the study of molecular systems. There are a number of excellent
simulation packages available, both free and commercial, which
incorporate many of these
methods.\cite{Brooks83,MacKerell98,Pearlman95,Berendsen95,Lindahl01,Smith96,Ponder87}
Because of our interest in rigid body dynamics, point multipoles, and
systems where the orientational degrees cannot be handled using the
common constraint procedures (like {\sc shake}), the majority of the
following work was performed using {\sc oopse}, the object-oriented
parallel simulation engine.\cite{Meineke05} The {\sc oopse} package
started out as a collection of separate programs written within our
group, and has developed into one of the few parallel molecular
dynamics packages capable of accurately integrating rigid bodies and
point multipoles.

In chapter \ref{chap:electrostatics}, we investigate correction
techniques for proper handling of long-ranged electrostatic
interactions. In particular we develop and evaluate some new pairwise
methods which we have incorporated into {\sc oopse}. These techniques
make an appearance in the later chapters, as they are applied to
specific systems of interest, showing how it they can improve the
quality of various molecular simulations.

In chapter \ref{chap:water}, we focus on simple water models,
specifically the single-point soft sticky dipole (SSD) family of water
models. These single-point models are more efficient than the common
multi-point partial charge models and, in many cases, better capture
the dynamic properties of water. We discuss improvements to these
models in regards to long-range electrostatic corrections and show
that these models can work well with the techniques discussed in
chapter \ref{chap:electrostatics}. By investigating and improving
simple water models, we are extending the limits of computational
efficiency for systems that heavily depend on water calculations.

In chapter \ref{chap:ice}, we study a unique polymorph of ice that we
discovered while performing water simulations with the fast simple
water models discussed in the previous chapter. This form of ice,
which we called ``imaginary ice'' (Ice-$i$), has a low-density
structure which is different from any known polymorph from either
experiment or other simulations. In this study, we perform a free
energy analysis and see that this structure is in fact the
thermodynamically preferred form of ice for both the single-point and
commonly used multi-point water models under the chosen simulation
conditions. We also consider electrostatic corrections, again
including the techniques discussed in chapter
\ref{chap:electrostatics}, to see how they influence the free energy
results. This work, to some degree, addresses the appropriateness of
using these simplistic water models outside of the conditions for
which they were developed.

In the final chapter we summarize the work presented previously. We
also close with some final comments before the supporting information
presented in the appendices.
Revision:	3019
Committed:	Fri Sep 22 13:45:24 2006 UTC (17 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	37451 byte(s)
Log Message:	more changes and an abstract included