trunk/chrisDissertation/Ice.tex

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

As discussed 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
property 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 and
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)
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 Ice-{\it i}
to indicate its origin in computational simulation. 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}A,B. These
tetramers form a crystal structure similar in appearance to a recent
two-dimensional surface tessellation simulated on silica.\cite{Yang04}
As expected in an ice crystal constructed of water tetramers, the
hydrogen bonds are not as linear as those observed 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-{\it i} and (B) Ice-$i^\prime$, the
elongated variant of Ice-{\it i}.  For Ice-{\it 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-{\it 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-{\it 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 work 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 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-{\it 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 ({\it NVT}) molecular dynamics calculations were
performed using the OOPSE molecular mechanics package.\cite{Meineke05}
The densities chosen for the simulations were taken from
isobaric-isothermal ({\it NPT}) simulations performed at 1 atm and
200~K. Each model (and each crystal structure) was allowed to relax for
300~ps in the {\it NPT} ensemble before averaging the density to obtain
the volumes for the {\it NVT} simulations.All molecules were treated
as rigid bodies, with orientational motion propagated using the
symplectic 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 over 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.  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
Fig. \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 has been shown to be reversible and provide
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
values 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 Particle-Mesh Ewald (PME) calculations for each
of the ice polymorphs.\cite{Ponder87} The calculated energy difference
in the presence and absence of PME 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 trucation 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
overwelmed 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}

The calculated free energies of proton-ordered variants of three low
density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-{\it 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-{\it i} is the most stable state for all of the
investigated water models.  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 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{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-{\it 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}

\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-{\it i} and the liquid state rather than ice I$_h$ and the liquid
state.  These results do not contradict other studies.  Studies of ice
I$_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-{\it i} were omitted, a $T_\textrm{m}$ value around 200~K
would be predicted from this work.  However, the $T_\textrm{m}$ from
Ice-{\it 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-{\it i} over all other polymorphs for this particular model under
these conditions.  While troubling, this behavior resulted in the
spontaneous crystallization of Ice-{\it 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}
\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-{\it 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
Fig. \ref{fig:incCutoff}, the free energies of the ice polymorphs with
water models lacking a long-range correction show significant cutoff
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-{\it 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 TIP3P was 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 particle mesh Ewald (PME).  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-{\it i}
polymorph.  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 preferred polymorph 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) & - \\
\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 ice polymorph Helmholtz free energies 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 it is not very close to Ice-$i$. The
SPC/E results really show the near isoenergetic behavior when using
the electrostatics 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, and the less
realistic ice B and Ice-$i^\prime$ 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.  However, the free
energy partially depends on simulation conditions (particularly on 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?  As
indicated above, 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 the identification of
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 Ice-$i$ 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, due to the stability of Ice-$i$ in the investigated
simulation conditions, the question arises as to possible experimental
observation of this polymorph.  The rather extensive past and current
experimental investigation of water in the low pressure regime makes
us hesitant to ascribe any relevance to this work outside of the
simulation community.  It is for this reason that we chose a name for
this polymorph which involves an imaginary quantity.  That said, there
are certain experimental conditions that would provide the most ideal
situation for possible observation. These include the negative
pressure or stretched solid regime, small clusters in vacuum
deposition environments, and in clathrate structures involving small
non-polar molecules.  For the purpose of comparison with experimental
results, we have calculated the oxygen-oxygen pair correlation
function, $g_\textrm{OO}(r)$, and the structure factor, $S(\vec{q})$
for the two Ice-{\it i} variants (along with example ice
I$_\textrm{h}$ and I$_\textrm{c}$ plots) at 77~K, and they are shown in
figures \ref{fig:gofr} and \ref{fig:sofq} respectively.  It is
interesting to note that the structure factors for Ice-$i^\prime$ and
Ice-I$_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 identify the new $i^\prime$ variant from
the I$_\textrm{c}$ polymorph.


\begin{figure}
\includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
\caption{Radial distribution functions of Ice-{\it 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-{\it 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. The labeled peaks compared favorably with ``spurious''
peaks observed in experimental studies of amorphous solid
water.\cite{Bizid87}}
\label{fig:sofq}
\end{figure}

Revision:	3019
Committed:	Fri Sep 22 13:45:24 2006 UTC (17 years, 11 months ago) by chrisfen
Content type:	application/x-tex
File size:	31086 byte(s)
Log Message:	more changes and an abstract included