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}

\section{\label{sec:MovingParticles}Propagating Particle Motion}

\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}, 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}.

\subsection{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}

\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{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}.

\section{Measuring Properties}

\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:	3001
Committed:	Thu Sep 7 20:42:27 2006 UTC (18 years, 7 months ago) by chrisfen
Content type:	application/x-tex
File size:	26493 byte(s)
Log Message:	Thar be changes in these documents