trunk/chrisDissertation/Introduction.tex

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

The following dissertation presents the primary aspects of the
research I have performed and been involved with over the last several
years. Rather than presenting the topics in a chronological fashion,
they were arranged to form a series where the later topics apply and
extend the findings of the former topics. This layout does lead to
occasional situations where knowledge gleaned from earlier chapters
(particularly chapter \ref{chap:electrostatics}) is not applied
outright in the later work; however, I feel that this organization is
more instructive and provides a more cohesive progression of research
efforts. 

This chapter gives a general overview of molecular simulations, with
particular emphasis on considerations that need to be made in order to
apply the technique properly. This leads quite naturally into chapter
\ref{chap:electrostatics}, where we investigate correction techniques
for proper handling of long-ranged electrostatic interactions. In
particular we develop and evaluate some new simple pairwise
methods. 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 depend heavily 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.

Finally, in chapter \ref{chap:conclusion}, we summarize the work
presented in the previous chapters and connect ideas together with
some final comments. The supporting information follows in the
appendix, and it gives a more detailed look at systems discussed in
chapter \ref{chap:electrostatics}.

\section{On Molecular Simulation}

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. 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 ($\mathbf{q}$) and velocity
($\dot{\mathbf{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(\mathbf{q},\dot{\mathbf{q}})dt = 0.
\end{equation}
The variation can be transferred to the variables that make up the
Lagrangian,
\begin{equation}
\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left(
        \frac{\partial L}{\partial \mathbf{q}_i}\delta \mathbf{q}_i
        + \frac{\partial L}{\partial \dot{\mathbf{q}}_i}\delta 
                \dot{\mathbf{q}}_i\right)dt = 0.
\end{equation}
Using the fact that $\dot{\mathbf{q}}$ is the derivative of
$\mathbf{q}$ with respect to time 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{\mathbf{q}}_i}
        - \frac{\partial L}{\partial \mathbf{q}_i}\right)
                \delta {\mathbf{q}}_i dt = 0,
\end{equation}
and since each variable is independent, we can separate the
contribution from each of the variables:
\begin{equation}
\frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i}
        - \frac{\partial L}{\partial \mathbf{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 can be calculated at each time
step with the 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 keep the system evolving in the region of phase-space
dictated by the ensemble and enjoy a greater degree of energy
conservation.\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
substantially larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope
{\it et al.}  developed a corrected version 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 forces are calculated at the new positions, 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_0, q_1, q_2,$ and $q_3$),
\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 method 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 our simulation code 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^3\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 {\sc dlm}
integration scheme.

\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 particles separated by less
than 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 particle within $R_\textrm{l}$ moves half the 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 the distance between a pair of particles fluctuates around
$R_\textrm{c}$, there will be sudden discontinuous jumps in the
potential (and forces) due to their inclusion and exclusion from the
interaction loop. 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 potential function
to eliminate these discontinuties, 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. 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}(r_{ij}) =
        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 inside $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 as well as in 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 interactions only from nearby
particles can lead to some issues, because particles at distances
greater than $R_\textrm{c}$ 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 interactions that extend 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 effect 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 properties can be
corrected by integration over the relevant function. Note that with
heterogeneous systems, this correction breaks 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 coordinates need not be
stored because they are easily calculated while determining particle
distances.

\section{Calculating Properties}

In order to use simulations to model experimental processes 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
mechanical 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 artificial 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. Since 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 for a
sufficiently chaotic system, and over a long enough period of time,
the time and ensemble averages are the same.

In addition to the average, the fluctuations of a particular property
can be determined via the standard deviation. Not only are
fluctuations useful for determining the spread of values around the
average and the error in the calculation of the value, they are also
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 other
thermodynamic properties, such as the constant pressure heat capacity
and the isothermal compressibility.

\section{OOPSE}

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 of freedom 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. This simulation package is
open-source software, available to anyone interested in performing
molecular dynamics simulations. More information about {\sc oopse} can
be found in reference \cite{Meineke05} or at the {\tt
http://oopse.org} website.


Revision:	3023
Committed:	Tue Sep 26 03:07:59 2006 UTC (17 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	39049 byte(s)
Log Message:	Refinements. Just a little proof checking left...