66 |
|
resulting in a pressure of approximately 1 atm at their respective |
67 |
|
temperatures. |
68 |
|
|
69 |
+ |
A single thermodynamic integration involves a sequence of simulations |
70 |
+ |
over which the system of interest is converted into a reference system |
71 |
+ |
for which the free energy is known. This transformation path is then |
72 |
+ |
integrated in order to determine the free energy difference between |
73 |
+ |
the two states: |
74 |
+ |
\begin{equation} |
75 |
+ |
\begin{center} |
76 |
+ |
\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda |
77 |
+ |
)}{\partial\lambda}\right\rangle_\lambda d\lambda, |
78 |
+ |
\end{center} |
79 |
+ |
\end{equation} |
80 |
+ |
where $V$ is the interaction potential and $\lambda$ is the |
81 |
+ |
transformation parameter. Simulations are distributed unevenly along |
82 |
+ |
this path in order to sufficiently sample the regions of greatest |
83 |
+ |
change in the potential. Typical integrations in this study consisted |
84 |
+ |
of $\sim$25 simulations ranging from 300 ps (for the unaltered system) |
85 |
+ |
to 75 ps (near the reference state) in length. |
86 |
+ |
|
87 |
|
For the thermodynamic integration of molecular crystals, the Einstein |
88 |
|
Crystal is chosen as the reference state that the system is converted |
89 |
|
to over the course of the simulation. In an Einstein Crystal, the |
259 |
|
risk of spontaneous crystallization. However, this risk changes when |
260 |
|
applying a longer cutoff. |
261 |
|
|
262 |
+ |
\begin{figure} |
263 |
+ |
\includegraphics[width=\linewidth]{cutoffChange.eps} |
264 |
+ |
\caption{Free energy as a function of cutoff radius for (A) SSD/E, (B) |
265 |
+ |
TIP3P, and (C) SSD/RF. Data points omitted include SSD/E: $I_c$ 12 |
266 |
+ |
\AA\, TIP3P: $I_c$ 12 \AA\ and B 12 \AA\, and SSD/RF: $I_c$ 9 |
267 |
+ |
\AA\. These crystals are unstable at 200 K and rapidly convert into a |
268 |
+ |
liquid. The connecting lines are qualitative visual aid.} |
269 |
+ |
\label{incCutoff} |
270 |
+ |
\end{figure} |
271 |
|
|
272 |
+ |
Increasing the cutoff radius in simulations of the more |
273 |
+ |
computationally efficient water models was done in order to evaluate |
274 |
+ |
the trend in free energy values when moving to systems that do not |
275 |
+ |
involve potential truncation. As seen in Fig. \ref{incCutoff}, the |
276 |
+ |
free energy of all the ice polymorphs show a substantial dependence on |
277 |
+ |
cutoff radius. In general, there is a narrowing of the free energy |
278 |
+ |
differences while moving to greater cutoff radius. This trend is much |
279 |
+ |
more subtle in the case of SSD/RF, indicating that the free energies |
280 |
+ |
calculated with a reaction field present provide a more accurate |
281 |
+ |
picture of the free energy landscape in the absence of potential |
282 |
+ |
truncation. |
283 |
|
|
284 |
+ |
To further study the changes resulting to the inclusion of a |
285 |
+ |
long-range interaction correction, the effect of an Ewald summation |
286 |
+ |
was estimated by applying the potential energy difference do to its |
287 |
+ |
inclusion in systems in the presence and absence of the |
288 |
+ |
correction. This was accomplished by calculation of the potential |
289 |
+ |
energy of identical crystals with and without PME using TINKER. The |
290 |
+ |
free energies for the investigated polymorphs using the TIP3P and |
291 |
+ |
SPC/E water models are shown in Table \ref{pmeShift}. TIP4P and TIP5P |
292 |
+ |
are not fully supported in TINKER, so the results for these models |
293 |
+ |
could not be estimated. The same trend pointed out through increase of |
294 |
+ |
cutoff radius is observed in these results. Ice-{\it i} is the |
295 |
+ |
preferred polymorph at ambient conditions for both the TIP3P and SPC/E |
296 |
+ |
water models; however, there is a narrowing of the free energy |
297 |
+ |
differences between the various solid forms. In the case of SPC/E this |
298 |
+ |
narrowing is significant enough that it becomes less clear cut that |
299 |
+ |
Ice-{\it i} is the most stable polymorph, and is possibly metastable |
300 |
+ |
with respect to ice B and possibly ice $I_c$. However, these results |
301 |
+ |
do not significantly alter the finding that the Ice-{\it i} polymorph |
302 |
+ |
is a stable crystal structure that should be considered when studying |
303 |
+ |
the phase behavior of water models. |
304 |
+ |
|
305 |
+ |
\begin{table*} |
306 |
+ |
\begin{minipage}{\linewidth} |
307 |
+ |
\renewcommand{\thefootnote}{\thempfootnote} |
308 |
+ |
\begin{center} |
309 |
+ |
\caption{The free energy of the studied ice polymorphs after applying |
310 |
+ |
the energy difference attributed to the inclusion of the PME |
311 |
+ |
long-range interaction correction. Units are kcal/mol.} |
312 |
+ |
\begin{tabular}{ l c c c c } |
313 |
+ |
\hline \\[-7mm] |
314 |
+ |
\ \ Water Model \ \ & \ \ \ \ \ $I_h$ \ \ & \ \ \ \ \ $I_c$ \ \ & \ \quad \ \ \ \ B \ \ & \ \ \ \ \ Ice-{\it i} \ \ \\ |
315 |
+ |
\hline \\[-3mm] |
316 |
+ |
\ \ TIP3P & \ \ -11.53 & \ \ -11.24 & \ \ -11.51 & \ \ -11.67\\ |
317 |
+ |
\ \ SPC/E & \ \ -12.77 & \ \ -12.92 & \ \ -12.96 & \ \ -13.02\\ |
318 |
+ |
\end{tabular} |
319 |
+ |
\label{pmeShift} |
320 |
+ |
\end{center} |
321 |
+ |
\end{minipage} |
322 |
+ |
\end{table*} |
323 |
+ |
|
324 |
|
\section{Conclusions} |
325 |
|
|
326 |
+ |
The free energy for proton ordered variants of hexagonal and cubic ice |
327 |
+ |
$I$, ice B, and recently discovered Ice-{\it i} where calculated under |
328 |
+ |
standard conditions for several common water models via thermodynamic |
329 |
+ |
integration. All the water models studied show Ice-{\it i} to be the |
330 |
+ |
minimum free energy crystal structure in the with a 9 \AA\ switching |
331 |
+ |
function cutoff. Calculated melting and boiling points show |
332 |
+ |
surprisingly good agreement with the experimental values; however, the |
333 |
+ |
solid phase at 1 atm is Ice-{\it i}, not ice $I_h$. The effect of |
334 |
+ |
interaction truncation was investigated through variation of the |
335 |
+ |
cutoff radius, use of a reaction field parameterized model, and |
336 |
+ |
estimation of the results in the presence of the Ewald summation |
337 |
+ |
correction. Interaction truncation has a significant effect on the |
338 |
+ |
computed free energy values, but Ice-{\it i} is still observed to be a |
339 |
+ |
relavent ice polymorph in simulation studies. |
340 |
+ |
|
341 |
|
\section{Acknowledgments} |
342 |
|
Support for this project was provided by the National Science |
343 |
|
Foundation under grant CHE-0134881. Computation time was provided by |
344 |
< |
the Notre Dame Bunch-of-Boxes (B.o.B) computer cluster under NSF grant |
345 |
< |
DMR-0079647. |
344 |
> |
the Notre Dame High Performance Computing Cluster and the Notre Dame |
345 |
> |
Bunch-of-Boxes (B.o.B) computer cluster (NSF grant DMR-0079647). |
346 |
|
|
347 |
|
\newpage |
348 |
|
|