| 1 |
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
| 2 |
\documentclass[11pt]{article} |
| 3 |
\usepackage{endfloat} |
| 4 |
\usepackage{amsmath} |
| 5 |
\usepackage{epsf} |
| 6 |
\usepackage{berkeley} |
| 7 |
\usepackage{setspace} |
| 8 |
\usepackage{tabularx} |
| 9 |
\usepackage{graphicx} |
| 10 |
\usepackage[ref]{overcite} |
| 11 |
\pagestyle{plain} |
| 12 |
\pagenumbering{arabic} |
| 13 |
\oddsidemargin 0.0cm \evensidemargin 0.0cm |
| 14 |
\topmargin -21pt \headsep 10pt |
| 15 |
\textheight 9.0in \textwidth 6.5in |
| 16 |
\brokenpenalty=10000 |
| 17 |
\renewcommand{\baselinestretch}{1.2} |
| 18 |
\renewcommand\citemid{\ } % no comma in optional reference note |
| 19 |
|
| 20 |
\begin{document} |
| 21 |
|
| 22 |
Canonical ensemble (NVT) molecular dynamics calculations were |
| 23 |
performed using the OOPSE molecular mechanics package.\cite{Meineke05} |
| 24 |
All molecules were treated as rigid bodies, with orientational motion |
| 25 |
propagated using the symplectic DLM integration method. Details about |
| 26 |
the implementation of this technique can be found in a recent |
| 27 |
publication.\cite{Dullweber1997} |
| 28 |
|
| 29 |
Thermodynamic integration is an established technique for |
| 30 |
determination of free energies of condensed phases of |
| 31 |
materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This |
| 32 |
method, implemented in the same manner illustrated by B\`{a}ez and |
| 33 |
Clancy, was utilized to calculate the free energy of several ice |
| 34 |
crystals at 200 K using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and |
| 35 |
SSD/E water models.\cite{Baez95a} Liquid state free energies at 300 |
| 36 |
and 400 K for all of these water models were also determined using |
| 37 |
this same technique, in order to determine melting points and to |
| 38 |
generate phase diagrams. System sizes were 648 or 1728 molecules for |
| 39 |
ice B, 1024 or 1280 molecules for ice $I_h$, 1000 molecules for ice |
| 40 |
$I_c$, and 1024 molecules for Ice-{\it i} and the liquid state |
| 41 |
simulations. The larger crystal sizes were necessary for simulations |
| 42 |
involving larger cutoff values. All simulations were carried out at |
| 43 |
densities which correspond to a pressure of approximately 1 atm at |
| 44 |
their respective temperatures. |
| 45 |
|
| 46 |
Thermodynamic integration involves a sequence of simulations during |
| 47 |
which the system of interest is converted into a reference system for |
| 48 |
which the free energy is known analytically. This transformation path |
| 49 |
is then integrated in order to determine the free energy difference |
| 50 |
between the two states: |
| 51 |
\begin{equation} |
| 52 |
\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda |
| 53 |
)}{\partial\lambda}\right\rangle_\lambda d\lambda, |
| 54 |
\end{equation} |
| 55 |
where $V$ is the interaction potential and $\lambda$ is the |
| 56 |
transformation parameter that scales the overall potential. |
| 57 |
Simulations are distributed strategically along this path in order to |
| 58 |
sufficiently sample the regions of greatest change in the potential. |
| 59 |
Typical integrations in this study consisted of $\sim$25 simulations |
| 60 |
ranging from 300 ps (for the unaltered system) to 75 ps (near the |
| 61 |
reference state) in length. |
| 62 |
|
| 63 |
For the thermodynamic integration of molecular crystals, the Einstein |
| 64 |
crystal was chosen as the reference system. In an Einstein crystal, |
| 65 |
the molecules are restrained at their ideal lattice locations and |
| 66 |
orientations. Using harmonic restraints, as applied by B\`{a}ez and |
| 67 |
Clancy, the total potential for this reference crystal |
| 68 |
($V_\mathrm{EC}$) is the sum of all the harmonic restraints, |
| 69 |
\begin{equation} |
| 70 |
V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} + |
| 71 |
\frac{K_\omega\omega^2}{2}, |
| 72 |
\end{equation} |
| 73 |
where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are |
| 74 |
the spring constants restraining translational motion and deflection |
| 75 |
of and rotation around the principle axis of the molecule |
| 76 |
respectively. These spring constants are typically calculated from |
| 77 |
the mean-square displacements of water molecules in an unrestrained |
| 78 |
ice crystal at 200 K. For these studies, $K_\mathrm{r} = 4.29$ kcal |
| 79 |
mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ = |
| 80 |
17.75$ kcal mol$^{-1}$. It is clear from Fig. \ref{waterSpring} that |
| 81 |
the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges |
| 82 |
from $-\pi$ to $\pi$. The partition function for a molecular crystal |
| 83 |
restrained in this fashion can be evaluated analytically, and the |
| 84 |
Helmholtz Free Energy ({\it A}) is given by |
| 85 |
\begin{eqnarray} |
| 86 |
A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left |
| 87 |
[\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right |
| 88 |
)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right |
| 89 |
)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right |
| 90 |
)^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi |
| 91 |
K_\omega K_\theta)^{\frac{1}{2}}}\exp\left |
| 92 |
(-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right |
| 93 |
)^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ], |
| 94 |
\label{ecFreeEnergy} |
| 95 |
\end{eqnarray} |
| 96 |
where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum |
| 97 |
potential energy of the ideal crystal.\cite{Baez95a} |
| 98 |
|
| 99 |
\begin{figure} |
| 100 |
\includegraphics[width=\linewidth]{rotSpring.eps} |
| 101 |
\caption{Possible orientational motions for a restrained molecule. |
| 102 |
$\theta$ angles correspond to displacement from the body-frame {\it |
| 103 |
z}-axis, while $\omega$ angles correspond to rotation about the |
| 104 |
body-frame {\it z}-axis. $K_\theta$ and $K_\omega$ are spring |
| 105 |
constants for the harmonic springs restraining motion in the $\theta$ |
| 106 |
and $\omega$ directions.} |
| 107 |
\label{waterSpring} |
| 108 |
\end{figure} |
| 109 |
|
| 110 |
In the case of molecular liquids, the ideal vapor is chosen as the |
| 111 |
target reference state. There are several examples of liquid state |
| 112 |
free energy calculations of water models present in the |
| 113 |
literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods |
| 114 |
typically differ in regard to the path taken for switching off the |
| 115 |
interaction potential to convert the system to an ideal gas of water |
| 116 |
molecules. In this study, we applied of one of the most convenient |
| 117 |
methods and integrated over the $\lambda^4$ path, where all |
| 118 |
interaction parameters are scaled equally by this transformation |
| 119 |
parameter. This method has been shown to be reversible and provide |
| 120 |
results in excellent agreement with other established |
| 121 |
methods.\cite{Baez95b} |
| 122 |
|
| 123 |
Charge, dipole, and Lennard-Jones interactions were modified by a |
| 124 |
cubic switching between 100\% and 85\% of the cutoff value (9 \AA). |
| 125 |
By applying this function, these interactions are smoothly truncated, |
| 126 |
thereby avoiding the poor energy conservation which results from |
| 127 |
harsher truncation schemes. The effect of a long-range correction was |
| 128 |
also investigated on select model systems in a variety of manners. |
| 129 |
For the SSD/RF model, a reaction field with a fixed dielectric |
| 130 |
constant of 80 was applied in all simulations.\cite{Onsager36} For a |
| 131 |
series of the least computationally expensive models (SSD/E, SSD/RF, |
| 132 |
TIP3P, and SPC/E), simulations were performed with longer cutoffs of |
| 133 |
10.5, 12, 13.5, and 15 \AA\ to compare with the 9 \AA\ cutoff results. |
| 134 |
Finally, the affects provided by an Ewald summation were estimated for |
| 135 |
TIP3P and SPC/E by performing single configuration calculations with |
| 136 |
Particle-Mesh Ewald (PME) in the TINKER molecular mechanics software |
| 137 |
package.\cite{Tinker} The calculated energy difference in the presence |
| 138 |
and absence of PME was applied to the previous results in order to |
| 139 |
predict changes to the free energy landscape. |
| 140 |
|
| 141 |
Additionally, $g_{OO}(r)$ and $S(\vec{q})$ plots were generated for |
| 142 |
the two Ice-{\it i} variants (along with example ice $I_h$ and $I_c$ |
| 143 |
plots) at 77K, and they are shown in figures \ref{fig:gofr} and |
| 144 |
\ref{fig:sofq}. The $S(\vec{q})$ is related to a three dimensional |
| 145 |
Fourier transform of the radial distribution function, which |
| 146 |
simplifies to the following expression: |
| 147 |
|
| 148 |
\begin{equation} |
| 149 |
S(q) = 1 + 4\pi\rho\int_{0}^{\infty} r^2 \frac{\sin kr}{kr}g_{OO}(r)dr, |
| 150 |
\label{sofqEq} |
| 151 |
\end{equation} |
| 152 |
|
| 153 |
where $\rho$ is the number density. To obtain a good estimation of |
| 154 |
$S(\vec{q})$, $g_{OO}(r)$ needs to extend to large $r$ values. Thus, |
| 155 |
simulations to obtain these plots were run on crystals eight times the |
| 156 |
size of those used in the thermodynamic integrations. |
| 157 |
|
| 158 |
\begin{figure} |
| 159 |
\includegraphics[width=\linewidth]{iceGofr.eps} |
| 160 |
\caption{Radial distribution functions of ice $I_h$, $I_c$, and |
| 161 |
Ice-{\it i} calculated from from simulations of the SSD/RF water model |
| 162 |
at 77 K. The Ice-{\it i} distribution function was obtained from |
| 163 |
simulations composed of TIP4P water.} |
| 164 |
\label{fig:gofr} |
| 165 |
\end{figure} |
| 166 |
|
| 167 |
\begin{figure} |
| 168 |
\includegraphics[width=\linewidth]{sofq.eps} |
| 169 |
\caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i}, |
| 170 |
and Ice-{\it i}$^\prime$ at 77 K. The raw structure factors have |
| 171 |
been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$ |
| 172 |
width) to compensate for the trunction effects in our finite size |
| 173 |
simulations.} |
| 174 |
\label{fig:sofq} |
| 175 |
\end{figure} |
| 176 |
|
| 177 |
|
| 178 |
\newpage |
| 179 |
|
| 180 |
\bibliographystyle{jcp} |
| 181 |
\bibliography{iceiSupplemental} |
| 182 |
|
| 183 |
|
| 184 |
\end{document} |
