trunk/chrisDissertation/Ice.tex

\chapter[PHASE BEHAVIOR OF WATER IN COMPUTER \\ SIMULATIONS]{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}

As mentioned in the previous chapter, water has proven to be a
challenging substance to depict in simulations, and a variety of
models have been developed to describe its behavior under varying
simulation
conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,vanderSpoel98,Urbic00,Mahoney00,Fennell04}
These models have been used to investigate important physical
phenomena like phase transitions and the hydrophobic
effect.\cite{Yamada02,Marrink94,Gallagher03} With the choice of models
available, it is only natural to compare them under interesting
thermodynamic conditions in an attempt to clarify the limitations of
each.\cite{Jorgensen83,Jorgensen98b,Baez94,Mahoney01} Two important
properties to quantify are the Gibbs and Helmholtz free energies,
particularly for the solid forms of water, as these predict the
thermodynamic stability of the various phases. Water has a
particularly rich phase diagram and takes on a number of different
stable crystalline structures as the temperature and pressure are
varied. This complexity makes it a challenging task to investigate the
entire free energy landscape.\cite{Sanz04} Ideally, research is
focused on the phases having the lowest free energy at a given state
point, because these phases will dictate the relevant transition
temperatures and pressures for the model.

The high-pressure phases of water (ice II-ice X as well as ice XII-ice XIV)
have been studied extensively both experimentally and
computationally. In this chapter, standard reference state methods
were applied in the {\it low} pressure regime to evaluate the free
energies for a few known crystalline water polymorphs that might be
stable at these pressures. This work is unique in the fact that one of
the crystal lattices was arrived at through crystallization of a
computationally efficient water model under constant pressure and
temperature conditions. 

While performing a series of melting simulations on an early iteration
of SSD/E, we observed several recrystallization events at a constant
pressure of 1 atm. After melting from ice I$_\textrm{h}$ at 235~K, two
of five systems recrystallized near 245~K. Crystallization events are
interesting in and of themselves;\cite{Matsumoto02,Yamada02} however,
the crystal structure extracted from these systems is different from
any previously observed ice polymorphs in experiment or
simulation.\cite{Fennell04} We have named this structure ``imaginary
ice'' (Ice-$i$) to indicate its origin in computer simulations. The
unit cell of Ice-$i$ and an axially elongated variant named
Ice-$i^\prime$ both consist of eight water molecules that stack in
rows of interlocking water tetramers as illustrated in figure
\ref{fig:iceiCell}. These tetramers form a crystal structure
similar in appearance to a recently simulated two-dimensional surface
tessellation of water on silica.\cite{Yang04} As expected in an ice
crystal constructed of water tetramers, the hydrogen bonds are not as
linear as those present in ice I$_\textrm{h}$; however, the
interlocking of these subunits appears to provide significant
stabilization to the overall crystal. The arrangement of these
tetramers results in open octagonal cavities that are typically
greater than 6.3~\AA\ in diameter (see figure
\ref{fig:protOrder}). This open structure leads to crystals that are
typically 0.07~g/cm$^3$ less dense than ice I$_\textrm{h}$.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/unitCell.pdf}
\caption{Unit cells for (A) Ice-$i$ and (B) Ice-$i^\prime$, the
elongated variant of Ice-$i$.  For Ice-$i$, the $a$ to $c$
relation is given by $a = 2.1214c$, while for Ice-$i^\prime$, $a =
1.7850c$.} 
\label{fig:iceiCell}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/orderedIcei.pdf}
\caption{Image of a proton ordered crystal of Ice-$i$ looking 
down the (001) crystal face. The rows of water tetramers surrounded by
octagonal pores leads to a crystal structure that is significantly
less dense than ice I$_\textrm{h}$.}
\label{fig:protOrder}
\end{figure}

Results from our initial studies indicated that Ice-$i$ is the
minimum energy crystal structure for the single point water models
investigated (for discussions on these single point dipole models, see
the previous chapter and related
articles\cite{Fennell04,Liu96,Bratko85}). These earlier results only
considered energetic stabilization and neglected entropic
contributions to the overall free energy. To address this issue, we
have calculated the absolute free energy of this crystal using
thermodynamic integration and compared it to the free energies of ice
I$_\textrm{c}$ and ice I$_\textrm{h}$ (the common low-density ice
polymorphs) and ice B (a higher density, but very stable crystal
structure observed by B\'{a}ez and Clancy in free energy studies of
SPC/E).\cite{Baez95b} This work includes results for the water model
from which Ice-$i$ was crystallized (SSD/E) in addition to several
common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
field parametrized single point dipole water model (SSD/RF). The
axially elongated variant, Ice-$i^\prime$, was used in calculations
involving SPC/E, TIP4P, and TIP5P. The square tetramer in Ice-$i$
distorts in Ice-$i^\prime$ to form a rhombus with alternating 85 and
95 degree angles. Under SPC/E, TIP4P, and TIP5P, this geometry is
better at forming favorable hydrogen bonds. The degree of rhomboid
distortion depends on the water model used but is significant enough
to split the peak in the radial distribution function which
corresponds to diagonal sites in the tetramers.

\section{Methods and Thermodynamic Integration}

Canonical ensemble ($NVT$) molecular dynamics calculations were
performed using the {\sc oopse} molecular mechanics
package.\cite{Meineke05} The densities chosen for the simulations were
taken from isobaric-isothermal ($NPT$) simulations performed at 1~atm
and 200~K. Each model (and each crystal structure) was allowed to
relax for 300~ps in the $NPT$ ensemble before averaging the density to
obtain the volumes for the $NVT$ simulations. All molecules were
treated as rigid bodies, with orientational motion propagated using
the symplectic {\sc dlm} integration method described in section
\ref{sec:IntroIntegrate}.


We used thermodynamic integration to calculate the Helmholtz free
energies ({\it A}) of the listed water models at various state
points. Thermodynamic integration is an established technique that has
been used extensively in the calculation of free energies for
condensed phases of
materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}.  This
method uses a sequence of simulations during which the system of
interest is converted into a reference system for which the free
energy is known analytically ($A_0$).  The difference in potential
energy between the reference system and the system of interest
($\Delta V$) is then integrated in order to determine the free energy
difference between the two states:
\begin{equation}
 A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda.
\end{equation}
Here, $\lambda$ is the parameter that governs the transformation
between the reference system and the system of interest.  For
crystalline phases, an harmonically-restrained (Einstein) crystal is
chosen as the reference state, while for liquid phases, the ideal gas
is taken as the reference state. Figure \ref{fig:integrationPath}
shows an example integration path for converting a crystalline system
to the Einstein crystal reference state.
\begin{figure}
\includegraphics[width=\linewidth]{./figures/integrationPath.pdf}
\caption{An example integration path to convert an unrestrained 
crystal ($\lambda = 1$) to the Einstein crystal reference state
($\lambda = 0$). Note the increase in samples at either end of the
path to improve the smoothness of the curve. For reversible processes,
conversion of the Einstein crystal back to the system of interest will
give an identical plot, thereby integrating to the same result.}
\label{fig:integrationPath}
\end{figure}

In an Einstein crystal, the molecules are restrained at their ideal
lattice locations and orientations. Using harmonic restraints, as
applied by B\'{a}ez and Clancy, the total potential for this reference
crystal ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
\begin{equation}
V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
\frac{K_\omega\omega^2}{2},
\end{equation}
where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
the spring constants restraining translational motion and deflection
of and rotation around the principle axis of the molecule respectively
(see figure \ref{fig:waterSpring}). These spring constants are
typically calculated from the mean-square displacements of water
molecules in an unrestrained ice crystal at 200~K.  For these studies,
$K_\mathrm{v} = 4.29$~kcal~mol$^{-1}$~\AA$^{-2}$, $K_\theta\ =
13.88$~kcal~mol$^{-1}$~rad$^{-2}$, and $K_\omega\ =
17.75$~kcal~mol$^{-1}$~rad$^{-2}$.  It is clear from figure
\ref{fig:waterSpring} that the values of $\theta$ range from $0$ to
$\pi$, while $\omega$ ranges from $-\pi$ to $\pi$.  The partition
function for a molecular crystal restrained in this fashion can be
evaluated analytically, and the Helmholtz Free Energy ({\it A}) is
given by
\begin{equation}
\begin{split}
A = E_m &- kT\ln\left(\frac{kT}{h\nu}\right)^3 \\
        &- kT\ln\left[\pi^\frac{1}{2}\left(
                \frac{8\pi^2I_\mathrm{A}kT}{h^2}\right)^\frac{1}{2}
                \left(\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right)^\frac{1}{2}
                \left(\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right)^\frac{1}{2}
                \right] \\ 
        &- kT\ln\left[\frac{kT}{2(\pi K_\omega K_\theta)^{\frac{1}{2}}}
                \exp\left(-\frac{kT}{2K_\theta}\right)
                \int_0^{\left(\frac{kT}{2K_\theta}\right)^\frac{1}{2}}
                \exp(t^2)\mathrm{d}t\right],
\end{split}
\label{eq:ecFreeEnergy}
\end{equation}
where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
potential energy of the ideal crystal.\cite{Baez95a} 

The choice of an Einstein crystal reference state is somewhat
arbitrary. Any ideal system for which the partition function is known
exactly could be used as a reference point, as long as the system does
not undergo a phase transition during the integration path between the
real and ideal systems.  Nada and van der Eerden have shown that the
use of different force constants in the Einstein crystal does not
affect the total free energy, and Gao {\it et al.} have shown that
free energies computed with the Debye crystal reference state differ
from the Einstein crystal by only a few tenths of a
kJ~mol$^{-1}$.\cite{Nada03,Gao00} These free energy differences can
lead to some uncertainty in the computed melting point of the solids.
\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/rotSpring.pdf}
\caption{Possible orientational motions for a restrained molecule. 
$\theta$ angles correspond to displacement from the body-frame {\it
z}-axis, while $\omega$ angles correspond to rotation about the
body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
constants for the harmonic springs restraining motion in the $\theta$
and $\omega$ directions.}
\label{fig:waterSpring}
\end{figure}

In the case of molecular liquids, the ideal vapor is chosen as the
target reference state.  There are several examples of liquid state
free energy calculations of water models present in the
literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
typically differ in regard to the path taken for switching off the
interaction potential to convert the system to an ideal gas of water
molecules.  In this study, we applied one of the most convenient
methods and integrated over the $\lambda^4$ path, where all
interaction parameters are scaled equally by this transformation
parameter.  This method is reversible and provides results in excellent
agreement with other established methods.\cite{Baez95b}

The Helmholtz free energy error was determined in the same manner in
both the solid and the liquid free energy calculations. At each point
along the integration path, we calculated the standard deviation of
the potential energy difference. Addition or subtraction of these
deviations to each of their respective points and integrating the
curve again provides the upper and lower bounds of the uncertainty in
the Helmholtz free energy.

Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
Lennard-Jones interactions were gradually reduced by a cubic switching
function.  By applying this function, these interactions are smoothly
truncated, thereby avoiding the poor energy conservation which results
from harsher truncation schemes.  The effect of a long-range
correction was also investigated on select model systems in a variety
of manners.  For the SSD/RF model, a reaction field with a fixed
dielectric constant of 80 was applied in all
simulations.\cite{Onsager36} For a series of the least computationally
expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were
performed with longer cutoffs of 10.5, 12, 13.5, and 15~\AA\ to
compare with the 9~\AA\ cutoff results.  Finally, the effects of using
the Ewald summation were estimated for TIP3P and SPC/E by performing
single configuration smooth particle-mesh Ewald ({\sc spme})
calculations for each of the ice polymorphs.\cite{Ponder87} The
calculated energy difference in the presence and absence of {\sc spme}
was applied to the previous results in order to predict changes to the
free energy landscape.

In addition to the above procedures, we also tested how the inclusion
of the Lennard-Jones long-range correction affects the free energy
results. The correction for the Lennard-Jones truncation was included
by integration of the equation discussed in section
\ref{sec:LJCorrections}. Rather than discuss its affect alongside the
free energy results, we will just mention that while the correction
does lower the free energy of the higher density states more than the
lower density states, the effect is so small that it is entirely
overwhelmed by the error in the free energy calculation. Since its
inclusion does not influence the results, the Lennard-Jones correction
was omitted from all the calculations below.

\section{Initial Free Energy Results}

\begin{table}
\centering
\caption{HELMHOLTZ FREE ENERGIES AND TRANSITION TEMPERATURES AT 1 
ATMOSPHERE FOR SEVERAL WATER MODELS}

\footnotesize
\begin{tabular}{lccccccc}
\toprule
\toprule
Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ & $T_\textrm{m}$ (*$T_\textrm{s}$) & $T_\textrm{b}$\\
\cmidrule(lr){2-6}
\cmidrule(l){7-8} 
& \multicolumn{5}{c}{(kcal mol$^{-1}$)} & \multicolumn{2}{c}{(K)}\\
\midrule
TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 266(7) & 337(4)\\
TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 262(6) & 354(4)\\
TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(7) & 357(4)\\
SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 299(6) & 396(4)\\
SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(4) & -\\
SSD/RF & -11.96(2) & -11.60(2) & -12.53(3) & -12.79(2) & - & 278(7) & 382(4)\\
\bottomrule
\end{tabular}
\label{tab:freeEnergy}
\end{table}
The calculated free energies of proton-ordered variants of three low
density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-$i$ or
Ice-$i^\prime$) and the stable higher density ice B are listed in
table \ref{tab:freeEnergy}.  Ice B was included because it has been
shown to be a minimum free energy structure for SPC/E at ambient
conditions.\cite{Baez95b} In addition to the free energies, the
relevant transition temperatures at standard pressure are also
displayed in table \ref{tab:freeEnergy}.  These free energy values
indicate that Ice-$i$ is the thermodynamically preferred state for
all of the investigated water models under the chosen simulation
conditions.  With the free energy at these state points, the
Gibbs-Helmholtz equation was used to project to other state points and
to build phase diagrams.  Figures \ref{fig:tp3PhaseDia} and
\ref{fig:ssdrfPhaseDia} are example phase diagrams built from the
results for the TIP3P and SSD/RF water models.  All other models have
similar structure, although the crossing points between the phases
move to different temperatures and pressures as indicated from the
transition temperatures in table
\ref{tab:freeEnergy}.  It is interesting to note that ice
I$_\textrm{h}$ (and ice I$_\textrm{c}$ for that matter) do not appear
in any of the phase diagrams for any of the models.  For purposes of
this study, ice B is representative of the dense ice polymorphs.  A
recent study by Sanz {\it et al.} provides details on the phase
diagrams for SPC/E and TIP4P at higher pressures than those studied
here.\cite{Sanz04}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/tp3PhaseDia.pdf}
\caption{Phase diagram for the TIP3P water model in the low pressure 
regime. The displayed $T_m$ and $T_b$ values are good predictions of
the experimental values; however, the solid phases shown are not the
experimentally observed forms. Both cubic and hexagonal ice I are
higher in energy and don't appear in the phase diagram.}
\label{fig:tp3PhaseDia}
\end{figure}
\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/ssdrfPhaseDia.pdf}
\caption{Phase diagram for the SSD/RF water model in the low pressure 
regime. Calculations producing these results were done under an
applied reaction field. It is interesting to note that this
computationally efficient model (over 3 times more efficient than
TIP3P) exhibits phase behavior similar to the less computationally
conservative charge based models.}
\label{fig:ssdrfPhaseDia}
\end{figure}

We note that all of the crystals investigated in this study are ideal
proton-ordered antiferroelectric structures. All of the structures
obey the Bernal-Fowler rules and should be able to form stable
proton-{\it disordered} crystals which have the traditional
$k_\textrm{B}$ln(3/2) residual entropy at 0~K.\cite{Bernal33,Pauling35}
Simulations of proton-disordered structures are relatively unstable
with all but the most expensive water models.\cite{Nada03} Our
simulations have therefore been performed with the ordered
antiferroelectric structures which do not require the residual entropy
term to be accounted for in the free energies. This may result in some
discrepancies when comparing our melting temperatures to the melting
temperatures that have been calculated via thermodynamic integrations
of the disordered structures.\cite{Sanz04}

Most of the water models have melting points that compare quite
favorably with the experimental value of 273~K.  The unfortunate
aspect of this result is that this phase change occurs between
Ice-$i$ and the liquid state rather than ice I$_\textrm{h}$ and the liquid
state.  These results do not contradict other studies.  Studies of ice
I$_\textrm{h}$ using TIP4P predict a $T_m$ ranging from 191 to 238~K
(differences being attributed to choice of interaction truncation and
different ordered and disordered molecular
arrangements).\cite{Nada03,Vlot99,Gao00,Sanz04} If the presence of ice
B and Ice-$i$ were omitted, a $T_\textrm{m}$ value around 200~K
would be predicted from this work.  However, the $T_\textrm{m}$ from
Ice-$i$ is calculated to be 262~K, indicating that these
simulation-based structures ought to be included in studies probing
phase transitions with this model.  Also of interest in these results
is that SSD/E does not exhibit a melting point at 1 atm but does
sublime at 355~K.  This is due to the significant stability of
Ice-$i$ over all other polymorphs for this particular model under
these conditions.  While troubling, this behavior resulted in the
spontaneous crystallization of Ice-$i$ which led us to investigate
this structure.  These observations provide a warning that simulations
of SSD/E as a ``liquid'' near 300~K are actually metastable and run
the risk of spontaneous crystallization.  However, when a longer
cutoff radius is used, SSD/E prefers the liquid state under standard
temperature and pressure.

\section{Effects of Potential Truncation}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
\caption{Free energy as a function of cutoff radius for SSD/E, TIP3P, 
SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
with an added Ewald correction term.  Error for the larger cutoff
points is equivalent to that observed at 9.0~\AA\ (see Table
\ref{tab:freeEnergy}). Data for ice I$_\textrm{c}$ with TIP3P using
both 12 and 13.5~\AA\ cutoffs were omitted because the crystal was
prone to distortion and melting at 200~K.  Ice-$i^\prime$ is the
form of Ice-$i$ used in the SPC/E simulations.}
\label{fig:incCutoff}
\end{figure}
For the more computationally efficient water models, we have also
investigated the effect of potential truncation on the computed free
energies as a function of the cutoff radius.  As seen in figure
\ref{fig:incCutoff}, the free energies of the ice polymorphs with
water models lacking a long-range correction show significant cutoff
radius dependence.  In general, there is a narrowing of the free
energy differences while moving to greater cutoff radii.  As the free
energies for the polymorphs converge, the stability advantage that
Ice-$i$ exhibits is reduced.  Adjacent to each of these plots are
results for systems with applied or estimated long-range corrections.
SSD/RF was parametrized for use with a reaction field, and the benefit
provided by this computationally inexpensive correction is apparent.
The free energies are largely independent of the size of the reaction
field cavity in this model, so small cutoff radii mimic bulk
calculations quite well under SSD/RF.
 
Although the point-charge models investigated here were parametrized
for use without the Ewald summation, we have estimated the effect of
this method for computing long-range electrostatics for both TIP3P and
SPC/E.  This was accomplished by calculating the potential energy of
identical crystals both with and without {\sc spme}.  Similar behavior
to that observed with reaction field is seen for both of these models.
The free energies show reduced dependence on cutoff radius and span a
narrower range for the various polymorphs.  Like the dipolar water
models, TIP3P displays a relatively constant preference for the
Ice-$i$ polymorph. The crystal preference is much more difficult to
determine for SPC/E.  Without a long-range correction, each of the
polymorphs studied assumes the role of the thermodynamically preferred
state under different cutoff radii.  The inclusion of the Ewald
correction flattens and narrows the gap in free energies such that the
polymorphs are isoenergetic within statistical uncertainty.  This
suggests that other conditions, such as the density in fixed-volume
simulations, can influence the polymorph expressed upon
crystallization.

\section{Expanded Results Using Damped Shifted Force Electrostatics}

In chapter \ref{chap:electrostatics}, we discussed in detail a
pairwise method for handling electrostatics (shifted force, {\sc sf})
that can be used as a simple and efficient replacement for the Ewald
summation. Answering the question of the free energies of these ice
polymorphs with varying water models would be an interesting
application of this technique. To this end, we set up thermodynamic
integrations of all of the previously discussed ice polymorphs using
the {\sc sf} technique with a cutoff radius of 12~\AA\ and an $\alpha$
of 0.2125~\AA . These calculations were performed on TIP5P-E and
TIP4P-Ew (variants of the root models optimized for the Ewald
summation) as well as SPC/E, SSD/RF, and TRED (see section
\ref{sec:tredWater}).

\begin{table}
\centering
\caption{HELMHOLTZ FREE ENERGIES OF ICE POLYMORPHS USING THE DAMPED
SHIFTED FORCE CORRECTION}
\begin{tabular}{ lccccc }
\toprule
\toprule
Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\
\cmidrule(lr){2-6}
& \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\
\midrule 
TIP5P-E & -11.98(4) & -11.96(4) & -11.87(3) & - & -11.95(3) \\
TIP4P-Ew & -13.11(3) & -13.09(3) & -12.97(3) & - & -12.98(3) \\
SPC/E & -12.99(3) & -13.00(3) & -13.03(3) & - & -12.99(3) \\
SSD/RF & -11.83(3) & -11.66(4) & -12.32(3) & -12.39(3) & - \\
TRED & -12.61(3) & -12.43(3) & -12.89(3) & -13.12(3) & - \\
\bottomrule
\end{tabular}
\label{tab:dampedFreeEnergy}
\end{table}
The results of these calculations in table \ref{tab:dampedFreeEnergy}
show similar behavior to the Ewald results in figure
\ref{fig:incCutoff}, at least for SSD/RF and SPC/E which are present
in both. The Helmholtz free energies of the ice polymorphs for SSD/RF
order in the same fashion; however, Ice-$i$ and ice B are quite a bit
closer in free energy (nearly isoenergetic). The free energy
differences between ice polymorphs for TRED water parallel SSD/RF,
with the exception that ice B is destabilized such that the free
energy is not very close to Ice-$i$. The SPC/E results really show the
near isoenergetic behavior when using the electrostatic
correction. Ice B has the lowest Helmholtz free energy; however, all
the polymorph results overlap within error.

The most interesting results from these calculations come from the
more expensive TIP4P-Ew and TIP5P-E results. Both of these models were
optimized for use with an electrostatic correction and are
geometrically arranged to mimic water following two different
ideas. In TIP5P-E, the primary location for the negative charge in the
molecule is assigned to the lone-pairs of the oxygen, while TIP4P-Ew
places the negative charge near the center-of-mass along the H-O-H
bisector. There is some debate as to which is the proper choice for
the negative charge location, and this has in part led to a six-site
water model that balances both of these options.\cite{Vega05,Nada03}
The limited results in table \ref{tab:dampedFreeEnergy} support the
results of Vega {\it et al.}, which indicate the TIP4P charge location
geometry is more physically valid.\cite{Vega05} With the TIP4P-Ew
water model, the experimentally observed polymorph (ice
I$_\textrm{h}$) is the preferred form with ice I$_\textrm{c}$ slightly
higher in energy, though overlapping within error. Additionally, the
less realistic ice B and Ice-$i^\prime$ structures are destabilized
relative to these polymorphs. TIP5P-E shows similar behavior to SPC/E,
where there is no real free energy distinction between the various
polymorphs, because many overlap within error. While ice B is close in
free energy to the other polymorphs, these results fail to support the
findings of other researchers indicating the preferred form of TIP5P
at 1~atm is a structure similar to ice
B.\cite{Yamada02,Vega05,Abascal05} It should be noted that we are
looking at TIP5P-E rather than TIP5P, and the differences in the
Lennard-Jones parameters could be a reason for this dissimilarity.
Overall, these results indicate that TIP4P-Ew is a better mimic of
real water than these other models when studying crystallization and
solid forms of water.

\section{Conclusions}

In this work, thermodynamic integration was used to determine the
absolute free energies of several ice polymorphs.  The new polymorph,
Ice-$i$ was observed to be the stable crystalline state for {\it all}
the water models when using a 9.0~\AA\ cutoff. Additional work showed
that the free energy depends in part on simulation conditions - in
particular, the choice of long-range correction method. Regardless,
Ice-$i$ was still observed to be a stable polymorph for all of the
studied water models.

So what is the preferred solid polymorph for simulated water?  The
answer appears to be dependent both on the conditions and the model
used.  In the case of short cutoffs without a long-range interaction
correction, Ice-$i$ and Ice-$i^\prime$ have the lowest free energy of
the studied polymorphs with all the models.  Ideally, crystallization
of each model under constant pressure conditions, as was done with
SSD/E, would aid in identifying their respective preferred structures.
This work, however, helps illustrate how studies involving one
specific model can lead to insight about important behavior of others.

We also note that none of the water models used in this study are
polarizable or flexible models. It is entirely possible that the
polarizability of real water makes the Ice-$i$ structure substantially
less stable than ice I$_\textrm{h}$. The dipole moment of the water
molecules increases as the system becomes more condensed, and the
increasing dipole moment should destabilize the tetramer structures in
Ice-$i$. Right now, using TIP4P-Ew with an electrostatic correction
gives the proper thermodynamically preferred state, and we recommend
this arrangement for study of crystallization processes if the
computational cost increase that comes with including polarizability
is an issue.

Finally, the stability of Ice-$i$ in these simulations raises the
possibility of experimental observation.  The extensive body of
research on water in the low pressure regime makes us hesitant to
ascribe any relevance to this work outside the simulation community.
It is for this reason that we chose a name for this polymorph
involving an imaginary quantity.  That said, there are certain
conditions that would be ideal for experimental observation of
Ice-$i$.  These include the negative pressure or stretched solid
regime, clusters deposited in vacuum environments, and clathrate
structures involving small non-polar molecules.  For the purpose of
comparison with future experimental results, we calculated the
oxygen-oxygen pair correlation function, $g_\textrm{OO}(r)$, and the
structure factor, $S(\vec{q})$ for the two Ice-$i$ variants (along
with ice I$_\textrm{h}$ and I$_\textrm{c}$) at 77~K (figures
\ref{fig:gofr} and \ref{fig:sofq}).  It is interesting to note that
the structure factors for Ice-$i^\prime$ and ice I$_\textrm{c}$ are
quite similar.  The primary differences are small peaks at 1.125,
2.29, and 2.53~\AA$^{-1}$, so particular attention to these regions
would be needed to distinguish Ice-$i^\prime$ from ice I$_\textrm{c}$.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
\caption{Radial distribution functions of Ice-$i$ and ice 
I$_\textrm{c}$ calculated from from simulations of the SSD/RF water
model at 77~K.}
\label{fig:gofr}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{./figures/sofq.pdf}
\caption{Predicted structure factors for Ice-$i$ and ice 
I$_\textrm{c}$ at 77~K.  The raw structure factors have been
convoluted with a Gaussian instrument function (0.075~\AA$^{-1}$
width) to compensate for the truncation effects in our finite size
simulations.}
\label{fig:sofq}
\end{figure}

Revision:	3045
Committed:	Fri Oct 13 21:03:14 2006 UTC (19 years ago) by chrisfen
Content type:	application/x-tex
File size:	30997 byte(s)
Log Message:	fixin figures and tables of contents