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 2979 by chrisfen, Tue Aug 29 00:40:05 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
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
13 > each of the models.\cite{Jorgensen83,Jorgensen98b,Baez94,Mahoney01}
14 > Two 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
17   particularly rich phase diagram and takes on a number of different and
# Line 23 | Line 24 | have been studied extensively both experimentally and
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 235K, two
38 > of five systems recrystallized near 245K. 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}$.
57  
# Line 68 | Line 74 | less dense than ice I$_\textrm{h}$.}
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
# 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}
# Line 104 | Line 110 | symplectic DLM integration method described in section
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} = 4.29$ kcal
163 > mol$^{-1}$ \AA$^{-2}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$ rad$^{-2}$,
164 > and $K_\omega\ = 17.75$ kcal mol$^{-1}$ rad$^{-2}$.  It is clear from
165 > Fig. \ref{fig:waterSpring} that the values of $\theta$ range from $0$ to
166 > $\pi$, while $\omega$ ranges from $-\pi$ to $\pi$.  The partition
167 > function for a molecular crystal restrained in this fashion can be
168 > evaluated analytically, and the Helmholtz Free Energy ({\it A}) is
169 > given by
170 > \begin{equation}
171 > \begin{split}
172 > A = E_m &- kT\ln\left(\frac{kT}{h\nu}\right)^3 \\
173 >        &- kT\ln\left[\pi^\frac{1}{2}\left(
174 >                \frac{8\pi^2I_\mathrm{A}kT}{h^2}\right)^\frac{1}{2}
175 >                \left(\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right)^\frac{1}{2}
176 >                \left(\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right)^\frac{1}{2}
177 >                \right] \\
178 >        &- kT\ln\left[\frac{kT}{2(\pi K_\omega K_\theta)^{\frac{1}{2}}}
179 >                \exp\left(-\frac{kT}{2K_\theta}\right)
180 >                \int_0^{\left(\frac{kT}{2K_\theta}\right)^\frac{1}{2}}
181 >                \exp(t^2)\mathrm{d}t\right],
182 > \end{split}
183   \label{eq:ecFreeEnergy}
184 < \end{eqnarray}
185 < where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$.\cite{Baez95a} In equation
186 < \ref{eq:ecFreeEnergy}, $K_\mathrm{v}$, $K_\mathrm{\theta}$, and
187 < $K_\mathrm{\omega}$ are the spring constants restraining translational
188 < motion and deflection of and rotation around the principle axis of the
189 < molecule respectively (Fig. \ref{waterSpring}), and $E_m$ is the
190 < minimum potential energy of the ideal crystal. In the case of
191 < molecular liquids, the ideal vapor is chosen as the target reference
192 < state.
193 <
184 > \end{equation}
185 > where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
186 > potential energy of the ideal crystal.\cite{Baez95a} The choice of an
187 > Einstein crystal reference stat is somewhat arbitrary. Any ideal
188 > system for which the partition function is known exactly could be used
189 > as a reference point as long as the system does not undergo a phase
190 > transition during the integration path between the real and ideal
191 > systems.  Nada and van der Eerden have shown that the use of different
192 > force constants in the Einstein crystal doesn not affect the total
193 > free energy, and Gao {\it et al.} have shown that free energies
194 > computed with the Debye crystal reference state differ from the
195 > Einstein crystal by only a few tenths of a kJ
196 > mol$^{-1}$.\cite{Nada03,Gao00} These free energy differences can lead
197 > to some uncertainty in the computed melting point of the solids.
198   \begin{figure}
199   \centering
200   \includegraphics[width=3.5in]{./figures/rotSpring.pdf}
# Line 164 | Line 204 | and $\omega$ directions.}
204   body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring
205   constants for the harmonic springs restraining motion in the $\theta$
206   and $\omega$ directions.}
207 < \label{waterSpring}
207 > \label{fig:waterSpring}
208   \end{figure}
209  
210 < Charge, dipole, and Lennard-Jones interactions were modified by a
211 < cubic switching between 100\% and 85\% of the cutoff value (9\AA). By
212 < applying this function, these interactions are smoothly truncated,
213 < thereby avoiding the poor energy conservation which results from
214 < harsher truncation schemes. The effect of a long-range correction was
215 < also investigated on select model systems in a variety of manners. For
216 < the SSD/RF model, a reaction field with a fixed dielectric constant of
217 < 80 was applied in all simulations.\cite{Onsager36} For a series of the
218 < least computationally expensive models (SSD/E, SSD/RF, and TIP3P),
219 < simulations were performed with longer cutoffs of 12 and 15\AA\ to
220 < compare with the 9\AA\ cutoff results. Finally, results from the use
221 < 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.
210 > In the case of molecular liquids, the ideal vapor is chosen as the
211 > target reference state.  There are several examples of liquid state
212 > free energy calculations of water models present in the
213 > literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
214 > typically differ in regard to the path taken for switching off the
215 > interaction potential to convert the system to an ideal gas of water
216 > molecules.  In this study, we applied one of the most convenient
217 > methods and integrated over the $\lambda^4$ path, where all
218 > interaction parameters are scaled equally by this transformation
219 > parameter.  This method has been shown to be reversible and provide
220 > results in excellent agreement with other established
221 > methods.\cite{Baez95b}
222  
223 < \section{Results and discussion}
223 > Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
224 > Lennard-Jones interactions were gradually reduced by a cubic switching
225 > function.  By applying this function, these interactions are smoothly
226 > truncated, thereby avoiding the poor energy conservation which results
227 > from harsher truncation schemes.  The effect of a long-range
228 > correction was also investigated on select model systems in a variety
229 > of manners.  For the SSD/RF model, a reaction field with a fixed
230 > dielectric constant of 80 was applied in all
231 > simulations.\cite{Onsager36} For a series of the least computationally
232 > expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were
233 > performed with longer cutoffs of 10.5, 12, 13.5, and 15\AA\ to
234 > compare with the 9\AA\ cutoff results.  Finally, the effects of using
235 > the Ewald summation were estimated for TIP3P and SPC/E by performing
236 > single configuration Particle-Mesh Ewald (PME) calculations for each
237 > of the ice polymorphs.\cite{Ponder87} The calculated energy difference
238 > in the presence and absence of PME was applied to the previous results
239 > in order to predict changes to the free energy landscape.
240  
241 < 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.
241 > \section{Initial Free Energy Results}
242  
243 + The calculated free energies of proton-ordered variants of three low
244 + density polymorphs (I$_\textrm{h}$, I$_\textrm{c}$, and Ice-{\it i} or
245 + Ice-$i^\prime$) and the stable higher density ice B are listed in
246 + table \ref{tab:freeEnergy}.  Ice B was included because it has been
247 + shown to be a minimum free energy structure for SPC/E at ambient
248 + conditions.\cite{Baez95b} In addition to the free energies, the
249 + relevant transition temperatures at standard pressure are also
250 + displayed in table \ref{tab:freeEnergy}.  These free energy values
251 + indicate that Ice-{\it i} is the most stable state for all of the
252 + investigated water models.  With the free energy at these state
253 + points, the Gibbs-Helmholtz equation was used to project to other
254 + state points and to build phase diagrams.  Figures
255 + \ref{fig:tp3PhaseDia} and \ref{fig:ssdrfPhaseDia} are example diagrams
256 + built from the results for the TIP3P and SSD/RF water models.  All
257 + other models have similar structure, although the crossing points
258 + between the phases move to different temperatures and pressures as
259 + indicated from the transition temperatures in table
260 + \ref{tab:freeEnergy}.  It is interesting to note that ice
261 + I$_\textrm{h}$ (and ice I$_\textrm{c}$ for that matter) do not appear
262 + in any of the phase diagrams for any of the models.  For purposes of
263 + this study, ice B is representative of the dense ice polymorphs.  A
264 + recent study by Sanz {\it et al.} provides details on the phase
265 + diagrams for SPC/E and TIP4P at higher pressures than those studied
266 + here.\cite{Sanz04}
267   \begin{table}
268   \centering
269 < \caption{HELMHOLTZ FREE ENERGIES FOR SEVERAL ICE POLYMORPHS WITH A
270 < VARIETY OF COMMON WATER MODELS AT 200 KELVIN AND 1 ATMOSPHERE}
271 < \begin{tabular}{ l  c  c  c  c }
269 > \caption{HELMHOLTZ FREE ENERGIES AND TRANSITION TEMPERATURES AT 1
270 > ATMOSPHERE FOR SEVERAL WATER MODELS}
271 >
272 > \footnotesize
273 > \begin{tabular}{lccccccc}
274   \toprule
275   \toprule
276 < Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or $i^\prime$) \\
277 < & kcal/mol & kcal/mol & kcal/mol & kcal/mol \\
276 > Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} & Ice-$i^\prime$ & $T_\textrm{m}$ (*$T_\textrm{s}$) & $T_\textrm{b}$\\
277 > \cmidrule(lr){2-6}
278 > \cmidrule(l){7-8}
279 > & \multicolumn{5}{c}{(kcal mol$^{-1}$)} & \multicolumn{2}{c}{(K)}\\
280   \midrule
281 < TIP3P & -11.41(4) & -11.23(6) & -11.82(5) & -12.30(5)\\
282 < TIP4P & -11.84(5) & -12.04(4) & -12.08(6) & -12.33(6)\\
283 < TIP5P & -11.85(5) & -11.86(4) & -11.96(4) & -12.29(4)\\
284 < SPC/E & -12.67(3) & -12.96(3) & -13.25(5) & -13.55(3)\\
285 < SSD/E & -11.27(3) & -11.19(8) & -12.09(4) & -12.54(4)\\
286 < SSD/RF & -11.51(4) & NA & -12.08(5) & -12.29(4)\\
281 > TIP3P & -11.41(2) & -11.23(3) & -11.82(3) & -12.30(3) & - & 269(7) & 357(4)\\
282 > TIP4P & -11.84(3) & -12.04(2) & -12.08(3) & - & -12.33(3) & 262(6) & 354(4)\\
283 > TIP5P & -11.85(3) & -11.86(2) & -11.96(2) & - & -12.29(2) & 266(7) & 337(4)\\
284 > SPC/E & -12.87(2) & -13.05(2) & -13.26(3) & - & -13.55(2) & 299(6) & 396(4)\\
285 > SSD/E & -11.27(2) & -11.19(4) & -12.09(2) & -12.54(2) & - & *355(4) & -\\
286 > SSD/RF & -11.96(2) & -11.60(2) & -12.53(3) & -12.79(2) & - & 278(7) & 382(4)\\
287   \bottomrule
288   \end{tabular}
289   \label{tab:freeEnergy}
290   \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}
291  
292   \begin{figure}
293   \centering
# Line 254 | Line 297 | higher in energy and don't appear in the phase diagram
297   the experimental values; however, the solid phases shown are not the
298   experimentally observed forms. Both cubic and hexagonal ice $I$ are
299   higher in energy and don't appear in the phase diagram.}
300 < \label{fig:tp3phasedia}
300 > \label{fig:tp3PhaseDia}
301   \end{figure}
302  
303   \begin{figure}
# Line 266 | Line 309 | conservative charge based models.}
309   computationally efficient model (over 3 times more efficient than
310   TIP3P) exhibits phase behavior similar to the less computationally
311   conservative charge based models.}
312 < \label{fig:ssdrfphasedia}
312 > \label{fig:ssdrfPhaseDia}
313   \end{figure}
314  
315 < \begin{table}
316 < \centering
317 < \caption{Melting ($T_m$), boiling ($T_b$), and sublimation ($T_s$)
318 < temperatures at 1 atm for several common water models compared with
319 < experiment. The $T_m$ and $T_s$ values from simulation correspond to a
320 < transition between Ice-{\it i} (or Ice-{\it i}$^\prime$) and the
321 < liquid or gas state.}
322 < \begin{tabular}{ l  c  c  c  c  c  c  c }
323 < \toprule
324 < \toprule
325 < Equilibria Point & TIP3P & TIP4P & TIP5P & SPC/E & SSD/E & SSD/RF & Exp.\\
326 < \midrule
327 < $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}
315 > We note that all of the crystals investigated in this study ar ideal
316 > proton-ordered antiferroelectric structures. All of the structures
317 > obey the Bernal-Fowler rules and should be able to form stable
318 > proton-{\it disordered} crystals which have the traditional
319 > $k_\textrm{B}$ln(3/2) residual entropy at 0K.\cite{Bernal33,Pauling35}
320 > Simulations of proton-disordered structures are relatively unstable
321 > with all but the most expensive water models.\cite{Nada03} Our
322 > simulations have therefore been performed with the ordered
323 > antiferroelectric structures which do not require the residual entropy
324 > term to be accounted for in the free energies. This may result in some
325 > discrepancies when comparing our melting temperatures to the melting
326 > temperatures that have been calculated via thermodynamic integrations
327 > of the disordered structures.\cite{Sanz04}
328  
329 < Table \ref{tab:meltandboil} lists the melting and boiling temperatures
330 < calculated from this work. Surprisingly, most of these models have
331 < melting points that compare quite favorably with experiment. The
332 < unfortunate aspect of this result is that this phase change occurs
333 < between Ice-{\it i} and the liquid state rather than ice
334 < 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
329 > Most of the water models have melting points that compare quite
330 > favorably with the experimental value of 273 K.  The unfortunate
331 > aspect of this result is that this phase change occurs between
332 > Ice-{\it i} and the liquid state rather than ice I$_h$ and the liquid
333 > state.  These results do not contradict other studies.  Studies of ice
334 > I$_h$ using TIP4P predict a $T_m$ ranging from 191 to 238 K
335   (differences being attributed to choice of interaction truncation and
336   different ordered and disordered molecular
337 < arrangements).\cite{Vlot99,Gao00,Sanz04} If the presence of ice B and
338 < Ice-{\it i} were omitted, a $T_m$ value around 210K would be predicted
339 < from this work. However, the $T_m$ from Ice-{\it i} is calculated at
340 < 265K, significantly higher in temperature than the previous
341 < studies. Also of interest in these results is that SSD/E does not
342 < exhibit a melting point at 1 atm, but it shows a sublimation point at
343 < 355K. This is due to the significant stability of Ice-{\it i} over all
344 < other polymorphs for this particular model under these
345 < conditions. While troubling, this behavior turned out to be
346 < advantageous in that it facilitated the spontaneous crystallization of
347 < Ice-{\it i}. These observations provide a warning that simulations of
348 < SSD/E as a ``liquid'' near 300K are actually metastable and run the
349 < risk of spontaneous crystallization. However, this risk changes when
350 < applying a longer cutoff.
337 > arrangements).\cite{Nada03,Vlot99,Gao00,Sanz04} If the presence of ice
338 > B and Ice-{\it i} were omitted, a $T_\textrm{m}$ value around 200 K
339 > would be predicted from this work.  However, the $T_\textrm{m}$ from
340 > Ice-{\it i} is calculated to be 262 K, indicating that these
341 > simulation based structures ought to be included in studies probing
342 > phase transitions with this model.  Also of interest in these results
343 > is that SSD/E does not exhibit a melting point at 1 atm but does
344 > sublime at 355 K.  This is due to the significant stability of
345 > Ice-{\it i} over all other polymorphs for this particular model under
346 > these conditions.  While troubling, this behavior resulted in the
347 > spontaneous crystallization of Ice-{\it i} which led us to investigate
348 > this structure.  These observations provide a warning that simulations
349 > of SSD/E as a ``liquid'' near 300 K are actually metastable and run
350 > the risk of spontaneous crystallization.  However, when a longer
351 > cutoff radius is used, SSD/E prefers the liquid state under standard
352 > temperature and pressure.
353  
354 + \section{Effects of Potential Trucation}
355 +
356   \begin{figure}
357   \includegraphics[width=\linewidth]{./figures/cutoffChange.pdf}
358 < \caption{Free energy as a function of cutoff radius for (A) SSD/E, (B)
359 < TIP3P, and (C) SSD/RF. Data points omitted include SSD/E:
360 < I$_\textrm{c}$ 12\AA\, TIP3P: I$_\textrm{c}$ 12\AA\ and B 12\AA\,
361 < and SSD/RF: I$_\textrm{c}$ 9\AA . These crystals are unstable at 200 K
362 < and rapidly convert into liquids. The connecting lines are qualitative
363 < visual aid.}
358 > \caption{Free energy as a function of cutoff radius for SSD/E, TIP3P,
359 > SPC/E, SSD/RF with a reaction field, and the TIP3P and SPC/E models
360 > with an added Ewald correction term.  Error for the larger cutoff
361 > points is equivalent to that observed at 9.0\AA\ (see Table
362 > \ref{tab:freeEnergy}). Data for ice I$_\textrm{c}$ with TIP3P using
363 > both 12 and 13.5\AA\ cutoffs were omitted because the crystal was
364 > prone to distortion and melting at 200K.  Ice-$i^\prime$ is the
365 > form of Ice-{\it i} used in the SPC/E simulations.}
366   \label{fig:incCutoff}
367   \end{figure}
368  
369 < Increasing the cutoff radius in simulations of the more
370 < computationally efficient water models was done in order to evaluate
371 < the trend in free energy values when moving to systems that do not
372 < involve potential truncation. As seen in figure \ref{fig:incCutoff},
373 < the free energy of all the ice polymorphs show a substantial
374 < dependence on cutoff radius. In general, there is a narrowing of the
375 < free energy differences while moving to greater cutoff
376 < radius. Interestingly, by increasing the cutoff radius, the free
377 < energy gap was narrowed enough in the SSD/E model that the liquid
378 < state is preferred under standard simulation conditions (298K and 1
379 < atm). Thus, it is recommended that simulations using this model choose
380 < interaction truncation radii greater than 9\AA\ . This narrowing
381 < trend is much more subtle in the case of SSD/RF, indicating that the
382 < free energies calculated with a reaction field present provide a more
383 < accurate picture of the free energy landscape in the absence of
384 < potential truncation.
369 > For the more computationally efficient water models, we have also
370 > investigated the effect of potential trunctaion on the computed free
371 > energies as a function of the cutoff radius.  As seen in
372 > Fig. \ref{fig:incCutoff}, the free energies of the ice polymorphs with
373 > water models lacking a long-range correction show significant cutoff
374 > dependence.  In general, there is a narrowing of the free energy
375 > differences while moving to greater cutoff radii.  As the free
376 > energies for the polymorphs converge, the stability advantage that
377 > Ice-{\it i} exhibits is reduced.  Adjacent to each of these plots are
378 > results for systems with applied or estimated long-range corrections.
379 > SSD/RF was parametrized for use with a reaction field, and the benefit
380 > provided by this computationally inexpensive correction is apparent.
381 > The free energies are largely independent of the size of the reaction
382 > field cavity in this model, so small cutoff radii mimic bulk
383 > calculations quite well under SSD/RF.
384 >
385 > Although TIP3P was paramaterized for use without the Ewald summation,
386 > we have estimated the effect of this method for computing long-range
387 > electrostatics for both TIP3P and SPC/E.  This was accomplished by
388 > calculating the potential energy of identical crystals both with and
389 > without particle mesh Ewald (PME).  Similar behavior to that observed
390 > with reaction field is seen for both of these models.  The free
391 > energies show reduced dependence on cutoff radius and span a narrower
392 > range for the various polymorphs.  Like the dipolar water models,
393 > TIP3P displays a relatively constant preference for the Ice-{\it i}
394 > polymorph.  Crystal preference is much more difficult to determine for
395 > SPC/E.  Without a long-range correction, each of the polymorphs
396 > studied assumes the role of the preferred polymorph under different
397 > cutoff radii.  The inclusion of the Ewald correction flattens and
398 > narrows the gap in free energies such that the polymorphs are
399 > isoenergetic within statistical uncertainty.  This suggests that other
400 > conditions, such as the density in fixed-volume simulations, can
401 > influence the polymorph expressed upon crystallization.
402  
403 < 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.
403 > \section{Expanded Results Using Damped Shifted Force Electrostatics}
404  
367 – \begin{table}
368 – \centering
369 – \caption{The free energy of the studied ice polymorphs after applying
370 – the energy difference attributed to the inclusion of the PME
371 – long-range interaction correction. Units are kcal/mol.}
372 – \begin{tabular}{ l  c  c  c  c }
373 – \toprule
374 – \toprule
375 – Water Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-{\it i} (or Ice-$i^\prime$) \\
376 – \midrule
377 – TIP3P  & -11.53(4) & -11.24(6) & -11.51(5) & -11.67(5)\\
378 – SPC/E  & -12.77(3) & -12.92(3) & -12.96(5) & -13.02(3)\\
379 – \bottomrule
380 – \end{tabular}
381 – \label{tab:pmeShift}
382 – \end{table}
405  
406   \section{Conclusions}
407  
408 < The free energy for proton ordered variants of hexagonal and cubic ice
409 < $I$, ice B, and recently discovered Ice-{\it i} were calculated under
410 < standard conditions for several common water models via thermodynamic
411 < integration. All the water models studied show Ice-{\it i} to be the
412 < minimum free energy crystal structure in the with a 9\AA\ switching
413 < function cutoff. Calculated melting and boiling points show
414 < surprisingly good agreement with the experimental values; however, the
415 < 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.
408 > In this work, thermodynamic integration was used to determine the
409 > absolute free energies of several ice polymorphs.  The new polymorph,
410 > Ice-{\it i} was observed to be the stable crystalline state for {\it
411 > all} the water models when using a 9.0\AA\ cutoff.  However, the free
412 > energy partially depends on simulation conditions (particularly on the
413 > choice of long range correction method). Regardless, Ice-{\it i} was
414 > still observered to be a stable polymorph for all of the studied water
415 > models.
416  
417 < Due to this relative stability of Ice-{\it i} in all manner of
418 < investigated simulation examples, the question arises as to possible
419 < experimental observation of this polymorph.  The rather extensive past
420 < and current experimental investigation of water in the low pressure
421 < regime makes us hesitant to ascribe any relevance of this work outside
422 < of the simulation community.  It is for this reason that we chose a
423 < name for this polymorph which involves an imaginary quantity.  That
424 < said, there are certain experimental conditions that would provide the
425 < most ideal situation for possible observation. These include the
426 < negative pressure or stretched solid regime, small clusters in vacuum
417 > So what is the preferred solid polymorph for simulated water?  As
418 > indicated above, the answer appears to be dependent both on the
419 > conditions and the model used.  In the case of short cutoffs without a
420 > long-range interaction correction, Ice-{\it i} and Ice-$i^\prime$ have
421 > the lowest free energy of the studied polymorphs with all the models.
422 > Ideally, crystallization of each model under constant pressure
423 > conditions, as was done with SSD/E, would aid in the identification of
424 > their respective preferred structures.  This work, however, helps
425 > illustrate how studies involving one specific model can lead to
426 > insight about important behavior of others.
427 >
428 > We also note that none of the water models used in this study are
429 > polarizable or flexible models.  It is entirely possible that the
430 > polarizability of real water makes Ice-{\it i} substantially less
431 > stable than ice I$_h$.  However, the calculations presented above seem
432 > interesting enough to communicate before the role of polarizability
433 > (or flexibility) has been thoroughly investigated.
434 >
435 > Finally, due to the stability of Ice-{\it i} in the investigated
436 > simulation conditions, the question arises as to possible experimental
437 > observation of this polymorph.  The rather extensive past and current
438 > experimental investigation of water in the low pressure regime makes
439 > us hesitant to ascribe any relevance to this work outside of the
440 > simulation community.  It is for this reason that we chose a name for
441 > this polymorph which involves an imaginary quantity.  That said, there
442 > are certain experimental conditions that would provide the most ideal
443 > situation for possible observation. These include the negative
444 > pressure or stretched solid regime, small clusters in vacuum
445   deposition environments, and in clathrate structures involving small
446 < non-polar molecules.  Figs. \ref{fig:gofr} and \ref{fig:sofq} contain
447 < our predictions for both the pair distribution function ($g_{OO}(r)$)
448 < and the structure factor ($S(\vec{q})$ for ice I$_\textrm{c}$ and for
449 < ice-{\it i} at a temperature of 77K.  In a quick comparison of the
450 < predicted S(q) for Ice-{\it i} and experimental studies of amorphous
451 < solid water, it is possible that some of the ``spurious'' peaks that
452 < could not be assigned in HDA could correspond to peaks labeled in this
453 < S(q).\cite{Bizid87} It should be noted that there is typically poor
454 < agreement on crystal densities between simulation and experiment, so
455 < such peak comparisons should be made with caution.  We will leave it
456 < to our experimental colleagues to determine whether this ice polymorph
425 < is named appropriately or if it should be promoted to Ice-0.
446 > non-polar molecules.  For the purpose of comparison with experimental
447 > results, we have calculated the oxygen-oxygen pair correlation
448 > function, $g_\textrm{OO}(r)$, and the structure factor, $S(\vec{q})$
449 > for the two Ice-{\it i} variants (along with example ice
450 > I$_\textrm{h}$ and I$_\textrm{c}$ plots) at 77K, and they are shown in
451 > figures \ref{fig:gofr} and \ref{fig:sofq} respectively.  It is
452 > interesting to note that the structure factors for Ice-$i^\prime$ and
453 > Ice-I$_c$ are quite similar.  The primary differences are small peaks
454 > at 1.125, 2.29, and 2.53\AA$^{-1}$, so particular attention to these
455 > regions would be needed to identify the new $i^\prime$ variant from
456 > the I$_\textrm{c}$ polymorph.
457  
458 +
459   \begin{figure}
460   \includegraphics[width=\linewidth]{./figures/iceGofr.pdf}
461   \caption{Radial distribution functions of Ice-{\it i} and ice

Diff Legend

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