ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/fennellDissertation/iceChapter.tex
(Generate patch)

Comparing trunk/fennellDissertation/iceChapter.tex (file contents):
Revision 2978 by chrisfen, Sun Aug 27 19:34:53 2006 UTC vs.
Revision 2986 by chrisfen, Wed Aug 30 22:14:37 2006 UTC

# Line 1 | Line 1
1   \chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}
2  
3 < Water has proven to be a challenging substance to depict in
4 < simulations, and a variety of models have been developed to describe
5 < its behavior under varying simulation
3 > As discussed in the previous chapter, water has proven to be a
4 > challenging substance to depict in simulations, and a variety of
5 > models have been developed to describe its behavior under varying
6 > simulation
7   conditions.\cite{Stillinger74,Rahman75,Berendsen81,Jorgensen83,Bratko85,Berendsen87,Caldwell95,Liu96,vanderSpoel98,Urbic00,Mahoney00,Fennell04}
8   These models have been used to investigate important physical
9   phenomena like phase transitions and the hydrophobic
10   effect.\cite{Yamada02,Marrink94,Gallagher03} With the choice of models
11 < available, it is only natural to compare the models under interesting
11 > available, it is only natural to compare them under interesting
12   thermodynamic conditions in an attempt to clarify the limitations of
13 < each of the models.\cite{Jorgensen83,Jorgensen98,Baez94,Mahoney01} Two
14 < important property to quantify are the Gibbs and Helmholtz free
15 < energies, particularly for the solid forms of water, as these predict
16 < the thermodynamic stability of the various phases. Water has a
13 > each.\cite{Jorgensen83,Jorgensen98b,Baez94,Mahoney01} Two important
14 > property to quantify are the Gibbs and Helmholtz free energies,
15 > particularly for the solid forms of water, as these predict the
16 > thermodynamic stability of the various phases. Water has a
17   particularly rich phase diagram and takes on a number of different and
18   stable crystalline structures as the temperature and pressure are
19   varied. This complexity makes it a challenging task to investigate the
# Line 23 | Line 24 | computatuionally. In this chapter, standard reference
24  
25   The high-pressure phases of water (ice II-ice X as well as ice XII)
26   have been studied extensively both experimentally and
27 < computatuionally. In this chapter, standard reference state methods
27 > computationally. In this chapter, standard reference state methods
28   were applied in the {\it low} pressure regime to evaluate the free
29   energies for a few known crystalline water polymorphs that might be
30   stable at these pressures. This work is unique in the fact that one of
31   the crystal lattices was arrived at through crystallization of a
32   computationally efficient water model under constant pressure and
33 < temperature conditions. Crystallization events are interesting in and
34 < of themselves;\cite{Matsumoto02,Yamada02} however, the crystal
35 < structure obtained in this case is different from any previously
36 < observed ice polymorphs in experiment or simulation.\cite{Fennell04}
37 < We have named this structure Ice-{\it i} to indicate its origin in
38 < computational simulation. The unit cell of Ice-$i$ and an axially
39 < elongated variant named Ice-$i^\prime$ both consist of eight water
40 < molecules that stack in rows of interlocking water tetramers as
41 < illustrated in figure \ref{fig:unitCell}A,B. These tetramers form a
42 < crystal structure similar in appearance to a recent two-dimensional
43 < surface tessellation simulated on silica.\cite{Yang04} As expected in
44 < an ice crystal constructed of water tetramers, the hydrogen bonds are
45 < not as linear as those observed in ice I$_\textrm{h}$; however, the
46 < interlocking of these subunits appears to provide significant
47 < stabilization to the overall crystal. The arrangement of these
48 < tetramers results in open octagonal cavities that are typically
49 < greater than 6.3\AA\ in diameter (see figure
33 > temperature conditions.
34 >
35 > While performing a series of melting simulations on an early iteration
36 > of SSD/E, we observed several recrystallization events at a constant
37 > pressure of 1 atm. After melting from ice I$_\textrm{h}$ at 235~K, two
38 > of five systems recrystallized near 245~K. Crystallization events are
39 > interesting in and of themselves;\cite{Matsumoto02,Yamada02} however,
40 > the crystal structure extracted from these systems is different from
41 > any previously observed ice polymorphs in experiment or
42 > simulation.\cite{Fennell04} We have named this structure Ice-{\it i}
43 > to indicate its origin in computational simulation. The unit cell of
44 > Ice-$i$ and an axially elongated variant named Ice-$i^\prime$ both
45 > consist of eight water molecules that stack in rows of interlocking
46 > water tetramers as illustrated in figure \ref{fig:iceiCell}A,B. These
47 > tetramers form a crystal structure similar in appearance to a recent
48 > two-dimensional surface tessellation simulated on silica.\cite{Yang04}
49 > As expected in an ice crystal constructed of water tetramers, the
50 > hydrogen bonds are not as linear as those observed in ice
51 > I$_\textrm{h}$; however, the interlocking of these subunits appears to
52 > provide significant stabilization to the overall crystal. The
53 > arrangement of these tetramers results in open octagonal cavities that
54 > are typically greater than 6.3~\AA\ in diameter (see figure
55   \ref{fig:protOrder}). This open structure leads to crystals that are
56 < typically 0.07 g/cm$^3$ less dense than ice I$_\textrm{h}$.
56 > typically 0.07~g/cm$^3$ less dense than ice I$_\textrm{h}$.
57  
58   \begin{figure}
59   \includegraphics[width=\linewidth]{./figures/unitCell.pdf}
# Line 68 | Line 74 | Results from our previous study indicated that Ice-{\i
74   \label{fig:protOrder}
75   \end{figure}
76  
77 < Results from our previous study indicated that Ice-{\it i} is the
77 > Results from our initial studies indicated that Ice-{\it i} is the
78   minimum energy crystal structure for the single point water models
79   investigated (for discussions on these single point dipole models, see
80   the previous work and related
81 < articles\cite{Fennell04,Liu96,Bratko85}). Our earlier results only
81 > articles\cite{Fennell04,Liu96,Bratko85}). These earlier results only
82   considered energetic stabilization and neglected entropic
83   contributions to the overall free energy. To address this issue, we
84   have calculated the absolute free energy of this crystal using
85   thermodynamic integration and compared to the free energies of ice
86   I$_\textrm{c}$ and ice I$_\textrm{h}$ (the common low-density ice
87   polymorphs) and ice B (a higher density, but very stable crystal
88 < structure observed by B\`{a}ez and Clancy in free energy studies of
88 > structure observed by B\'{a}ez and Clancy in free energy studies of
89   SPC/E).\cite{Baez95b} This work includes results for the water model
90   from which Ice-{\it i} was crystallized (SSD/E) in addition to several
91   common water models (TIP3P, TIP4P, TIP5P, and SPC/E) and a reaction
# Line 93 | Line 99 | to diagonal sites in the tetramers.
99   to split the peak in the radial distribution function which corresponds
100   to diagonal sites in the tetramers.
101  
102 < \section{Methods}
102 > \section{Methods and Thermodynamic Integration}
103  
104   Canonical ensemble ({\it NVT}) molecular dynamics calculations were
105   performed using the OOPSE molecular mechanics package.\cite{Meineke05}
106   The densities chosen for the simulations were taken from
107   isobaric-isothermal ({\it NPT}) simulations performed at 1 atm and
108 < 200K. Each model (and each crystal structure) was allowed to relax for
109 < 300ps in the {\it NPT} ensemble before averaging the density to obtain
108 > 200~K. Each model (and each crystal structure) was allowed to relax for
109 > 300~ps in the {\it NPT} ensemble before averaging the density to obtain
110   the volumes for the {\it NVT} simulations.All molecules were treated
111   as rigid bodies, with orientational motion propagated using the
112   symplectic DLM integration method described in section
113 < \ref{sec:IntroIntegration}.
113 > \ref{sec:IntroIntegrate}.
114  
115 +
116   We used thermodynamic integration to calculate the Helmholtz free
117   energies ({\it A}) of the listed water models at various state
118   points. Thermodynamic integration is an established technique that has
119   been used extensively in the calculation of free energies for
120   condensed phases of
121 < materials.\cite{Frenkel84,Hermans88,Meijer90,Baez95,Vlot99} This
122 < method uses a sequence of simulations during which the system of
121 > materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}.  This
122 > method uses a sequence of simulations over which the system of
123   interest is converted into a reference system for which the free
124 < energy is known analytically ($A_0$). This transformation path is then
125 < integrated in order to determine the free energy difference between
126 < the two states:
124 > energy is known analytically ($A_0$).  The difference in potential
125 > energy between the reference system and the system of interest
126 > ($\Delta V$) is then integrated in order to determine the free energy
127 > difference between the two states:
128   \begin{equation}
129 < \Delta A = \int_0^1\left\langle\frac{\partial V(\lambda
122 < )}{\partial\lambda}\right\rangle_\lambda d\lambda,
129 > A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda.
130   \end{equation}
131 < where $V$ is the interaction potential and $\lambda$ is the
132 < transformation parameter that scales the overall
133 < potential. Simulations are distributed unevenly along this path in
134 < order to sufficiently sample the regions of greatest change in the
135 < potential. Typical integrations in this study consisted of $\sim$25
136 < simulations ranging from 300ps (for the unaltered system) to 75ps
137 < (near the reference state) in length.
131 > Here, $\lambda$ is the parameter that governs the transformation
132 > between the reference system and the system of interest.  For
133 > crystalline phases, an harmonically-restrained (Einstein) crystal is
134 > chosen as the reference state, while for liquid phases, the ideal gas
135 > is taken as the reference state. Figure \ref{fig:integrationPath}
136 > shows an example integration path for converting a crystalline system
137 > to the Einstein crystal reference state.
138 > \begin{figure}
139 > \includegraphics[width=\linewidth]{./figures/integrationPath.pdf}
140 > \caption{An example integration path to convert an unrestrained
141 > crystal ($\lambda = 1$) to the Einstein crystal reference state
142 > ($\lambda = 0$). Note the increase in samples at either end of the
143 > path to improve the smoothness of the curve. For reversible processes,
144 > conversion of the Einstein crystal back to the system of interest will
145 > give an identical plot, thereby integrating to the same result.}
146 > \label{fig:integrationPath}
147 > \end{figure}
148  
149 < For the thermodynamic integration of molecular crystals, the Einstein
150 < crystal was chosen as the reference state. In an Einstein crystal, the
151 < molecules are harmonically restrained at their ideal lattice locations
152 < and orientations. The partition function for a molecular crystal
153 < restrained in this fashion can be evaluated analytically, and the
154 < Helmholtz Free Energy ({\it A}) is given by
155 < \begin{eqnarray}
156 < A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
157 < [\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
158 < )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
159 < )^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
160 < )^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
161 < K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
162 < (-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
163 < )^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
149 > In an Einstein crystal, the molecules are restrained at their ideal
150 > lattice locations and orientations. Using harmonic restraints, as
151 > applied by B\'{a}ez and Clancy, the total potential for this reference
152 > crystal ($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
153 > \begin{equation}
154 > V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
155 > \frac{K_\omega\omega^2}{2},
156 > \end{equation}
157 > where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
158 > the spring constants restraining translational motion and deflection
159 > of and rotation around the principle axis of the molecule
160 > respectively.  These spring constants are typically calculated from
161 > the mean-square displacements of water molecules in an unrestrained
162 > ice crystal at 200~K.  For these studies, $K_\mathrm{v} =
163 > 4.29$~kcal~mol$^{-1}$~\AA$^{-2}$, $K_\theta\ =
164 > 13.88$~kcal~mol$^{-1}$~rad$^{-2}$, and $K_\omega\ =
165 > 17.75$~kcal~mol$^{-1}$~rad$^{-2}$.  It is clear from
166 > Fig. \ref{fig:waterSpring} that the values of $\theta$ range from $0$
167 > to $\pi$, while $\omega$ ranges from $-\pi$ to $\pi$.  The partition
168 > function for a molecular crystal restrained in this fashion can be
169 > evaluated analytically, and the Helmholtz Free Energy ({\it A}) is
170 > given by
171 > \begin{equation}
172 > \begin{split}
173 > A = E_m &- kT\ln\left(\frac{kT}{h\nu}\right)^3 \\
174 >        &- kT\ln\left[\pi^\frac{1}{2}\left(
175 >                \frac{8\pi^2I_\mathrm{A}kT}{h^2}\right)^\frac{1}{2}
176 >                \left(\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right)^\frac{1}{2}
177 >                \left(\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right)^\frac{1}{2}
178 >                \right] \\
179 >        &- kT\ln\left[\frac{kT}{2(\pi K_\omega K_\theta)^{\frac{1}{2}}}
180 >                \exp\left(-\frac{kT}{2K_\theta}\right)
181 >                \int_0^{\left(\frac{kT}{2K_\theta}\right)^\frac{1}{2}}
182 >                \exp(t^2)\mathrm{d}t\right],
183 > \end{split}
184   \label{eq:ecFreeEnergy}
185 < \end{eqnarray}
186 < where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
187 < \ref{eq:ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
188 < $K_\mathrm{\omega}$ are the spring constants restraining translational
189 < motion and deflection of and rotation around the principle axis of the
190 < molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
191 < minimum potential energy of the ideal crystal. In the case of
192 < molecular liquids, the ideal vapor is chosen as the target reference
193 < state.
194 <
185 > \end{equation}
186 > where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
187 > potential energy of the ideal crystal.\cite{Baez95a} The choice of an
188 > Einstein crystal reference state is somewhat arbitrary. Any ideal
189 > system for which the partition function is known exactly could be used
190 > as a reference point as long as the system does not undergo a phase
191 > transition during the integration path between the real and ideal
192 > systems.  Nada and van der Eerden have shown that the use of different
193 > force constants in the Einstein crystal does not affect the total
194 > free energy, and Gao {\it et al.} have shown that free energies
195 > computed with the Debye crystal reference state differ from the
196 > Einstein crystal by only a few tenths of a
197 > kJ~mol$^{-1}$.\cite{Nada03,Gao00} These free energy differences can
198 > lead to some uncertainty in the computed melting point of the solids.
199   \begin{figure}
200   \centering
201   \includegraphics[width=3.5in]{./figures/rotSpring.pdf}
# Line 164 | Line 205 | and $\omega$ directions.}
205   body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
206   constants for the harmonic springs restraining motion in the $\theta$
207   and $\omega$ directions.}
208 < \label{waterSpring}
208 > \label{fig:waterSpring}
209   \end{figure}
210  
211 < Charge, dipole, and Lennard-Jones interactions were modified by a
212 < cubic switching between 100\% and 85\% of the cutoff value (9\AA). By
213 < applying this function, these interactions are smoothly truncated,
214 < thereby avoiding the poor energy conservation which results from
215 < harsher truncation schemes. The effect of a long-range correction was
216 < also investigated on select model systems in a variety of manners. For
217 < the SSD/RF model, a reaction field with a fixed dielectric constant of
218 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
219 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
220 < simulations were performed with longer cutoffs of 12 and 15\AA\ to
221 < compare with the 9\AA\ cutoff results. Finally, results from the use
222 < of an Ewald summation were estimated for TIP3P and SPC/E by performing
182 < calculations with Particle-Mesh Ewald (PME) in the TINKER molecular
183 < mechanics software package.\cite{Tinker} The calculated energy
184 < difference in the presence and absence of PME was applied to the
185 < previous results in order to predict changes to the free energy
186 < landscape.
211 > In the case of molecular liquids, the ideal vapor is chosen as the
212 > target reference state.  There are several examples of liquid state
213 > free energy calculations of water models present in the
214 > literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
215 > typically differ in regard to the path taken for switching off the
216 > interaction potential to convert the system to an ideal gas of water
217 > molecules.  In this study, we applied one of the most convenient
218 > methods and integrated over the $\lambda^4$ path, where all
219 > interaction parameters are scaled equally by this transformation
220 > parameter.  This method has been shown to be reversible and provide
221 > results in excellent agreement with other established
222 > methods.\cite{Baez95b}
223  
224 < \section{Results and discussion}
224 > Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
225 > Lennard-Jones interactions were gradually reduced by a cubic switching
226 > function.  By applying this function, these interactions are smoothly
227 > truncated, thereby avoiding the poor energy conservation which results
228 > from harsher truncation schemes.  The effect of a long-range
229 > correction was also investigated on select model systems in a variety
230 > of manners.  For the SSD/RF model, a reaction field with a fixed
231 > dielectric constant of 80 was applied in all
232 > simulations.\cite{Onsager36} For a series of the least computationally
233 > expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were
234 > performed with longer cutoffs of 10.5, 12, 13.5, and 15~\AA\ to
235 > compare with the 9~\AA\ cutoff results.  Finally, the effects of using
236 > the Ewald summation were estimated for TIP3P and SPC/E by performing
237 > single configuration Particle-Mesh Ewald (PME) calculations for each
238 > of the ice polymorphs.\cite{Ponder87} The calculated energy difference
239 > in the presence and absence of PME was applied to the previous results
240 > in order to predict changes to the free energy landscape.
241  
242 < The free energy of proton ordered Ice-{\it i} was calculated and
191 < compared with the free energies of proton ordered variants of the
192 < experimentally observed low-density ice polymorphs, I$_\textrm{h}$ and
193 < I$_\textrm{c}$, as well as the higher density ice B, observed by
194 < B\`{a}ez and Clancy and thought to be the minimum free energy
195 < structure for the SPC/E model at ambient conditions.\cite{Baez95b} Ice
196 < XI, the experimentally-observed proton-ordered variant of ice
197 < I$_\textrm{h}$, was investigated initially, but was found to be not as
198 < stable as proton disordered or antiferroelectric variants of ice
199 < I$_\textrm{h}$. The proton ordered variant of ice I$_\textrm{h}$ used
200 < here is a simple antiferroelectric version that has an 8 molecule unit
201 < cell.\cite{Davidson84} The crystals contained 648 or 1728 molecules
202 < for ice B, 1024 or 1280 molecules for ice I$_\textrm{h}$, 1000
203 < molecules for ice I$_\textrm{c}$, or 1024 molecules for Ice-{\it
204 < i}. The larger crystal sizes were necessary for simulations involving
205 < larger cutoff values.
242 > \section{Initial Free Energy Results}
243  
244 + The calculated free energies of proton-ordered variants of three low
245 + density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-{\it i} or
246 + Ice-$i^\prime$) and the stable higher density ice B are listed in
247 + table \ref{tab:freeEnergy}.  Ice B was included because it has been
248 + shown to be a minimum free energy structure for SPC/E at ambient
249 + conditions.\cite{Baez95b} In addition to the free energies, the
250 + relevant transition temperatures at standard pressure are also
251 + displayed in table \ref{tab:freeEnergy}.  These free energy values
252 + indicate that Ice-{\it i} is the most stable state for all of the
253 + investigated water models.  With the free energy at these state
254 + points, the Gibbs-Helmholtz equation was used to project to other
255 + state points and to build phase diagrams.  Figures
256 + \ref{fig:tp3PhaseDia} and \ref{fig:ssdrfPhaseDia} are example diagrams
257 + built from the results for the TIP3P and SSD/RF water models.  All
258 + other models have similar structure, although the crossing points
259 + between the phases move to different temperatures and pressures as
260 + indicated from the transition temperatures in table
261 + \ref{tab:freeEnergy}.  It is interesting to note that ice
262 + I$_\textrm{h}$ (and ice I$_\textrm{c}$ for that matter) do not appear
263 + in any of the phase diagrams for any of the models.  For purposes of
264 + this study, ice B is representative of the dense ice polymorphs.  A
265 + recent study by Sanz {\it et al.} provides details on the phase
266 + diagrams for SPC/E and TIP4P at higher pressures than those studied
267 + here.\cite{Sanz04}
268   \begin{table}
269   \centering
270 < \caption{HELMHOLTZ FREE ENERGIES FOR SEVERAL ICE POLYMORPHS WITH A
271 < VARIETY OF COMMON WATER MODELS AT 200 KELVIN AND 1 ATMOSPHERE}
272 < \begin{tabular}{ l  c  c  c  c }
270 > \caption{HELMHOLTZ FREE ENERGIES AND TRANSITION TEMPERATURES AT 1
271 > ATMOSPHERE FOR SEVERAL WATER MODELS}
272 >
273 > \footnotesize
274 > \begin{tabular}{lccccccc}
275   \toprule
276   \toprule
277 < Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or $i^\prime$) \\
278 < & kcal/mol & kcal/mol & kcal/mol & kcal/mol \\
277 > Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} & Ice-$i^\prime$ & $T_\textrm{m}$ (*$T_\textrm{s}$) & $T_\textrm{b}$\\
278 > \cmidrule(lr){2-6}
279 > \cmidrule(l){7-8}
280 > & \multicolumn{5}{c}{(kcal mol$^{-1}$)} & \multicolumn{2}{c}{(K)}\\
281   \midrule
282 < TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\
283 < TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\
284 < TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\
285 < SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\
286 < SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\
287 < SSD/RF & -11.51(4) & NA & -12.08(5) & -12.29(4)\\
282 > TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 266(7) & 337(4)\\
283 > TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 262(6) & 354(4)\\
284 > TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(7) & 357(4)\\
285 > SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 299(6) & 396(4)\\
286 > SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(4) & -\\
287 > SSD/RF & -11.96(2) & -11.60(2) & -12.53(3) & -12.79(2) & - & 278(7) & 382(4)\\
288   \bottomrule
289   \end{tabular}
290   \label{tab:freeEnergy}
291   \end{table}
227 –
228 – Table \ref{tab:freeEnergy} shows the results of the free energy
229 – calculations with a cutoff radius of 9\AA. It should be noted that the
230 – ice I$_\textrm{c}$ crystal polymorph is not stable at 200K and 1 atm
231 – with the SSD/RF water model, hense omitted results for that cell. The
232 – free energy values displayed in this table, it is clear that Ice-{\it
233 – i} (or Ice-$i^\prime$ for TIP4P, TIP5P, and SPC/E) is the most stable
234 – state for all of the common water models studied.
235 –
236 – With the free energy at these state points, the Gibbs-Helmholtz
237 – equation was used to project to other state points and to build phase
238 – diagrams.  Figures \ref{fig:tp3phasedia} and \ref{fig:ssdrfphasedia}
239 – are example diagrams built from the free energy results. All other
240 – models have similar structure, although the crossing points between
241 – the phases exist at slightly different temperatures and pressures. It
242 – is interesting to note that ice I does not exist in either cubic or
243 – hexagonal form in any of the phase diagrams for any of the models. For
244 – purposes of this study, ice B is representative of the dense ice
245 – polymorphs. A recent study by Sanz {\it et al.} goes into detail on
246 – the phase diagrams for SPC/E and TIP4P in the high pressure
247 – regime.\cite{Sanz04}
292  
293   \begin{figure}
294   \centering
# Line 254 | Line 298 | higher in energy and don't appear in the phase diagram
298   the experimental values; however, the solid phases shown are not the
299   experimentally observed forms. Both cubic and hexagonal ice $I$ are
300   higher in energy and don't appear in the phase diagram.}
301 < \label{fig:tp3phasedia}
301 > \label{fig:tp3PhaseDia}
302   \end{figure}
303  
304   \begin{figure}
# Line 266 | Line 310 | conservative charge based models.}
310   computationally efficient model (over 3 times more efficient than
311   TIP3P) exhibits phase behavior similar to the less computationally
312   conservative charge based models.}
313 < \label{fig:ssdrfphasedia}
313 > \label{fig:ssdrfPhaseDia}
314   \end{figure}
315  
316 < \begin{table}
317 < \centering
318 < \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
319 < temperatures at 1 atm for several common water models compared with
320 < experiment. The $T_m$ and $T_s$ values from simulation correspond to a
321 < transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
322 < liquid or gas state.}
323 < \begin{tabular}{ l  c  c  c  c  c  c  c }
324 < \toprule
325 < \toprule
326 < Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
327 < \midrule
328 < $T_m$ (K)  & 269(8) & 266(10) & 271(7) & 296(5) & - & 278(7) & 273\\
285 < $T_b$ (K)  & 357(2) & 354(3) & 337(3) & 396(4) & - & 348(3) & 373\\
286 < $T_s$ (K)  & - & - & - & - & 355(3) & - & -\\
287 < \bottomrule
288 < \end{tabular}
289 < \label{tab:meltandboil}
290 < \end{table}
316 > We note that all of the crystals investigated in this study are ideal
317 > proton-ordered antiferroelectric structures. All of the structures
318 > obey the Bernal-Fowler rules and should be able to form stable
319 > proton-{\it disordered} crystals which have the traditional
320 > $k_\textrm{B}$ln(3/2) residual entropy at 0~K.\cite{Bernal33,Pauling35}
321 > Simulations of proton-disordered structures are relatively unstable
322 > with all but the most expensive water models.\cite{Nada03} Our
323 > simulations have therefore been performed with the ordered
324 > antiferroelectric structures which do not require the residual entropy
325 > term to be accounted for in the free energies. This may result in some
326 > discrepancies when comparing our melting temperatures to the melting
327 > temperatures that have been calculated via thermodynamic integrations
328 > of the disordered structures.\cite{Sanz04}
329  
330 < Table \ref{tab:meltandboil} lists the melting and boiling temperatures
331 < calculated from this work. Surprisingly, most of these models have
332 < melting points that compare quite favorably with experiment. The
333 < unfortunate aspect of this result is that this phase change occurs
334 < between Ice-{\it i} and the liquid state rather than ice
335 < I$_\textrm{h}$ and the liquid state. These results are actually not
298 < contrary to previous studies in the literature. Earlier free energy
299 < studies of ice I using TIP4P predict a $T_m$ ranging from 214 to 238K
330 > Most of the water models have melting points that compare quite
331 > favorably with the experimental value of 273~K.  The unfortunate
332 > aspect of this result is that this phase change occurs between
333 > Ice-{\it i} and the liquid state rather than ice I$_h$ and the liquid
334 > state.  These results do not contradict other studies.  Studies of ice
335 > I$_h$ using TIP4P predict a $T_m$ ranging from 191 to 238~K
336   (differences being attributed to choice of interaction truncation and
337   different ordered and disordered molecular
338 < arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
339 < Ice-{\it i} were omitted, a $T_m$ value around 210K would be predicted
340 < from this work. However, the $T_m$ from Ice-{\it i} is calculated at
341 < 265K, significantly higher in temperature than the previous
342 < studies. Also of interest in these results is that SSD/E does not
343 < exhibit a melting point at 1 atm, but it shows a sublimation point at
344 < 355K. This is due to the significant stability of Ice-{\it i} over all
345 < other polymorphs for this particular model under these
346 < conditions. While troubling, this behavior turned out to be
347 < advantageous in that it facilitated the spontaneous crystallization of
348 < Ice-{\it i}. These observations provide a warning that simulations of
349 < SSD/E as a ``liquid'' near 300K are actually metastable and run the
350 < risk of spontaneous crystallization. However, this risk changes when
351 < applying a longer cutoff.
338 > arrangements).\cite{Nada03,Vlot99,Gao00,Sanz04} If the presence of ice
339 > B and Ice-{\it i} were omitted, a $T_\textrm{m}$ value around 200~K
340 > would be predicted from this work.  However, the $T_\textrm{m}$ from
341 > Ice-{\it i} is calculated to be 262~K, indicating that these
342 > simulation based structures ought to be included in studies probing
343 > phase transitions with this model.  Also of interest in these results
344 > is that SSD/E does not exhibit a melting point at 1 atm but does
345 > sublime at 355~K.  This is due to the significant stability of
346 > Ice-{\it i} over all other polymorphs for this particular model under
347 > these conditions.  While troubling, this behavior resulted in the
348 > spontaneous crystallization of Ice-{\it i} which led us to investigate
349 > this structure.  These observations provide a warning that simulations
350 > of SSD/E as a ``liquid'' near 300~K are actually metastable and run
351 > the risk of spontaneous crystallization.  However, when a longer
352 > cutoff radius is used, SSD/E prefers the liquid state under standard
353 > temperature and pressure.
354  
355 + \section{Effects of Potential Truncation}
356 +
357   \begin{figure}
358   \includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
359 < \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
360 < TIP3P, and (C) SSD/RF. Data points omitted include SSD/E:
361 < I$_\textrm{c}$ 12\AA\, TIP3P: I$_\textrm{c}$ 12\AA\ and B 12\AA\,
362 < and SSD/RF: I$_\textrm{c}$ 9\AA . These crystals are unstable at 200 K
363 < and rapidly convert into liquids. The connecting lines are qualitative
364 < visual aid.}
359 > \caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
360 > SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
361 > with an added Ewald correction term.  Error for the larger cutoff
362 > points is equivalent to that observed at 9.0~\AA\ (see Table
363 > \ref{tab:freeEnergy}). Data for ice I$_\textrm{c}$ with TIP3P using
364 > both 12 and 13.5~\AA\ cutoffs were omitted because the crystal was
365 > prone to distortion and melting at 200~K.  Ice-$i^\prime$ is the
366 > form of Ice-{\it i} used in the SPC/E simulations.}
367   \label{fig:incCutoff}
368   \end{figure}
369  
370 < Increasing the cutoff radius in simulations of the more
371 < computationally efficient water models was done in order to evaluate
372 < the trend in free energy values when moving to systems that do not
373 < involve potential truncation. As seen in figure \ref{fig:incCutoff},
374 < the free energy of all the ice polymorphs show a substantial
375 < dependence on cutoff radius. In general, there is a narrowing of the
376 < free energy differences while moving to greater cutoff
377 < radius. Interestingly, by increasing the cutoff radius, the free
378 < energy gap was narrowed enough in the SSD/E model that the liquid
379 < state is preferred under standard simulation conditions (298K and 1
380 < atm). Thus, it is recommended that simulations using this model choose
381 < interaction truncation radii greater than 9\AA\ . This narrowing
382 < trend is much more subtle in the case of SSD/RF, indicating that the
383 < free energies calculated with a reaction field present provide a more
384 < accurate picture of the free energy landscape in the absence of
385 < potential truncation.
370 > For the more computationally efficient water models, we have also
371 > investigated the effect of potential truncation on the computed free
372 > energies as a function of the cutoff radius.  As seen in
373 > Fig. \ref{fig:incCutoff}, the free energies of the ice polymorphs with
374 > water models lacking a long-range correction show significant cutoff
375 > dependence.  In general, there is a narrowing of the free energy
376 > differences while moving to greater cutoff radii.  As the free
377 > energies for the polymorphs converge, the stability advantage that
378 > Ice-{\it i} exhibits is reduced.  Adjacent to each of these plots are
379 > results for systems with applied or estimated long-range corrections.
380 > SSD/RF was parametrized for use with a reaction field, and the benefit
381 > provided by this computationally inexpensive correction is apparent.
382 > The free energies are largely independent of the size of the reaction
383 > field cavity in this model, so small cutoff radii mimic bulk
384 > calculations quite well under SSD/RF.
385 >
386 > Although TIP3P was parametrized for use without the Ewald summation,
387 > we have estimated the effect of this method for computing long-range
388 > electrostatics for both TIP3P and SPC/E.  This was accomplished by
389 > calculating the potential energy of identical crystals both with and
390 > without particle mesh Ewald (PME).  Similar behavior to that observed
391 > with reaction field is seen for both of these models.  The free
392 > energies show reduced dependence on cutoff radius and span a narrower
393 > range for the various polymorphs.  Like the dipolar water models,
394 > TIP3P displays a relatively constant preference for the Ice-{\it i}
395 > polymorph.  Crystal preference is much more difficult to determine for
396 > SPC/E.  Without a long-range correction, each of the polymorphs
397 > studied assumes the role of the preferred polymorph under different
398 > cutoff radii.  The inclusion of the Ewald correction flattens and
399 > narrows the gap in free energies such that the polymorphs are
400 > isoenergetic within statistical uncertainty.  This suggests that other
401 > conditions, such as the density in fixed-volume simulations, can
402 > influence the polymorph expressed upon crystallization.
403  
404 < To further study the changes resulting to the inclusion of a
346 < long-range interaction correction, the effect of an Ewald summation
347 < was estimated by applying the potential energy difference do to its
348 < inclusion in systems in the presence and absence of the
349 < correction. This was accomplished by calculation of the potential
350 < energy of identical crystals with and without PME using TINKER. The
351 < free energies for the investigated polymorphs using the TIP3P and
352 < SPC/E water models are shown in Table \ref{tab:pmeShift}. TIP4P and
353 < TIP5P are not fully supported in TINKER, so the results for these
354 < models could not be estimated. The same trend pointed out through
355 < increase of cutoff radius is observed in these PME results. Ice-{\it
356 < i} is the preferred polymorph at ambient conditions for both the TIP3P
357 < and SPC/E water models; however, there is a narrowing of the free
358 < energy differences between the various solid forms. In the case of
359 < SPC/E this narrowing is significant enough that it becomes less clear
360 < that Ice-{\it i} is the most stable polymorph, and is possibly
361 < metastable with respect to ice B and possibly ice
362 < I$_\textrm{c}$. However, these results do not significantly alter the
363 < finding that the Ice-{\it i} polymorph is a stable crystal structure
364 < that should be considered when studying the phase behavior of water
365 < models.
404 > \section{Expanded Results Using Damped Shifted Force Electrostatics}
405  
406 + In chapter \ref{chap:electrostatics}, we discussed in detail a
407 + pairwise method for handling electrostatics (shifted force, {\sc sf})
408 + that can be used as a simple and efficient replacement for the Ewald
409 + summation. Answering the question of the free energies of these ice
410 + polymorphs with varying water models would be an interesting
411 + application of this technique. To this end, we set up thermodynamic
412 + integrations of all of the previously discussed ice polymorphs using
413 + the {\sc sf} technique with a cutoff radius of 12~\AA\ and an $\alpha$
414 + of 0.2125~\AA . These calculations were performed on TIP5P-E and
415 + TIP4P-Ew (variants of the root models optimized for the Ewald
416 + summation) as well as SPC/E, SSD/RF, and TRED (see section
417 + \ref{sec:tredWater}).
418 +
419   \begin{table}
420   \centering
421 < \caption{The free energy of the studied ice polymorphs after applying
422 < the energy difference attributed to the inclusion of the PME
423 < long-range interaction correction. Units are kcal/mol.}
372 < \begin{tabular}{ l  c  c  c  c }
421 > \caption{HELMHOLTZ FREE ENERGIES OF ICE POLYMORPHS USING THE DAMPED
422 > SHIFTED FORCE CORRECTION}
423 > \begin{tabular}{ lccccc }
424   \toprule
425   \toprule
426 < Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or Ice-$i^\prime$) \\
427 < \midrule
428 < TIP3P  & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\
429 < SPC/E  & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\
430 < \bottomrule
426 > Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\
427 > \cmidrule(lr){2-6}
428 > & \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\
429 > \midrule
430 > TIP5P-E & -10.76(4) & -10.72(4) & & - & -10.68(4) \\
431 > TIP4P-Ew & & -11.77(3) & & - & -11.60(3) \\
432 > SPC/E & -12.98(3) & -11.60(3) & & - & -12.93(3) \\
433 > SSD/RF & -11.81(4) & -11.65(3) & & -12.41(4) & - \\
434 > TRED & -12.58(3) & -12.44(3) & & -13.09(4) & - \\
435   \end{tabular}
436 < \label{tab:pmeShift}
436 > \label{tab:dampedFreeEnergy}
437   \end{table}
438  
439 +
440   \section{Conclusions}
441  
442 < The free energy for proton ordered variants of hexagonal and cubic ice
443 < $I$, ice B, and recently discovered Ice-{\it i} were calculated under
444 < standard conditions for several common water models via thermodynamic
445 < integration. All the water models studied show Ice-{\it i} to be the
446 < minimum free energy crystal structure in the with a 9\AA\ switching
447 < function cutoff. Calculated melting and boiling points show
448 < surprisingly good agreement with the experimental values; however, the
449 < solid phase at 1 atm is Ice-{\it i}, not ice I$_\textrm{h}$. The
394 < effect of interaction truncation was investigated through variation of
395 < the cutoff radius, use of a reaction field parameterized model, and
396 < estimation of the results in the presence of the Ewald
397 < summation. Interaction truncation has a significant effect on the
398 < computed free energy values, and may significantly alter the free
399 < energy landscape for the more complex multipoint water models. Despite
400 < these effects, these results show Ice-{\it i} to be an important ice
401 < polymorph that should be considered in simulation studies.
442 > In this work, thermodynamic integration was used to determine the
443 > absolute free energies of several ice polymorphs.  The new polymorph,
444 > Ice-{\it i} was observed to be the stable crystalline state for {\it
445 > all} the water models when using a 9.0~\AA\ cutoff.  However, the free
446 > energy partially depends on simulation conditions (particularly on the
447 > choice of long range correction method). Regardless, Ice-{\it i} was
448 > still observed to be a stable polymorph for all of the studied water
449 > models.
450  
451 < Due to this relative stability of Ice-{\it i} in all manner of
452 < investigated simulation examples, the question arises as to possible
453 < experimental observation of this polymorph.  The rather extensive past
454 < and current experimental investigation of water in the low pressure
455 < regime makes us hesitant to ascribe any relevance of this work outside
456 < of the simulation community.  It is for this reason that we chose a
457 < name for this polymorph which involves an imaginary quantity.  That
458 < said, there are certain experimental conditions that would provide the
459 < most ideal situation for possible observation. These include the
460 < negative pressure or stretched solid regime, small clusters in vacuum
451 > So what is the preferred solid polymorph for simulated water?  As
452 > indicated above, the answer appears to be dependent both on the
453 > conditions and the model used.  In the case of short cutoffs without a
454 > long-range interaction correction, Ice-{\it i} and Ice-$i^\prime$ have
455 > the lowest free energy of the studied polymorphs with all the models.
456 > Ideally, crystallization of each model under constant pressure
457 > conditions, as was done with SSD/E, would aid in the identification of
458 > their respective preferred structures.  This work, however, helps
459 > illustrate how studies involving one specific model can lead to
460 > insight about important behavior of others.
461 >
462 > We also note that none of the water models used in this study are
463 > polarizable or flexible models.  It is entirely possible that the
464 > polarizability of real water makes Ice-{\it i} substantially less
465 > stable than ice I$_h$.  However, the calculations presented above seem
466 > interesting enough to communicate before the role of polarizability
467 > (or flexibility) has been thoroughly investigated.
468 >
469 > Finally, due to the stability of Ice-{\it i} in the investigated
470 > simulation conditions, the question arises as to possible experimental
471 > observation of this polymorph.  The rather extensive past and current
472 > experimental investigation of water in the low pressure regime makes
473 > us hesitant to ascribe any relevance to this work outside of the
474 > simulation community.  It is for this reason that we chose a name for
475 > this polymorph which involves an imaginary quantity.  That said, there
476 > are certain experimental conditions that would provide the most ideal
477 > situation for possible observation. These include the negative
478 > pressure or stretched solid regime, small clusters in vacuum
479   deposition environments, and in clathrate structures involving small
480 < non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
481 < our predictions for both the pair distribution function ($g_{OO}(r)$)
482 < and the structure factor ($S(\vec{q})$ for ice I$_\textrm{c}$ and for
483 < ice-{\it i} at a temperature of 77K.  In a quick comparison of the
484 < predicted S(q) for Ice-{\it i} and experimental studies of amorphous
485 < solid water, it is possible that some of the ``spurious'' peaks that
486 < could not be assigned in HDA could correspond to peaks labeled in this
487 < S(q).\cite{Bizid87} It should be noted that there is typically poor
488 < agreement on crystal densities between simulation and experiment, so
489 < such peak comparisons should be made with caution.  We will leave it
490 < to our experimental colleagues to determine whether this ice polymorph
425 < is named appropriately or if it should be promoted to Ice-0.
480 > non-polar molecules.  For the purpose of comparison with experimental
481 > results, we have calculated the oxygen-oxygen pair correlation
482 > function, $g_\textrm{OO}(r)$, and the structure factor, $S(\vec{q})$
483 > for the two Ice-{\it i} variants (along with example ice
484 > I$_\textrm{h}$ and I$_\textrm{c}$ plots) at 77~K, and they are shown in
485 > figures \ref{fig:gofr} and \ref{fig:sofq} respectively.  It is
486 > interesting to note that the structure factors for Ice-$i^\prime$ and
487 > Ice-I$_c$ are quite similar.  The primary differences are small peaks
488 > at 1.125, 2.29, and 2.53~\AA$^{-1}$, so particular attention to these
489 > regions would be needed to identify the new $i^\prime$ variant from
490 > the I$_\textrm{c}$ polymorph.
491  
492 +
493   \begin{figure}
494   \includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
495   \caption{Radial distribution functions of Ice-{\it i} and ice
496   I$_\textrm{c}$ calculated from from simulations of the SSD/RF water
497 < model at 77 K.}
497 > model at 77~K.}
498   \label{fig:gofr}
499   \end{figure}
500  
501   \begin{figure}
502   \includegraphics[width=\linewidth]{./figures/sofq.pdf}
503   \caption{Predicted structure factors for Ice-{\it i} and ice
504 < I$_\textrm{c}$ at 77 K.  The raw structure factors have been
505 < convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
506 < width) to compensate for the trunction effects in our finite size
504 > I$_\textrm{c}$ at 77~K.  The raw structure factors have been
505 > convoluted with a gaussian instrument function (0.075~\AA$^{-1}$
506 > width) to compensate for the truncation effects in our finite size
507   simulations. The labeled peaks compared favorably with ``spurious''
508   peaks observed in experimental studies of amorphous solid
509   water.\cite{Bizid87}}

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines