trunk/scienceIcei/iceiSupplemental.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
\documentclass[11pt]{article}
%\usepackage{endfloat}
\usepackage{amsmath}
\usepackage{epsf}
\usepackage{berkeley}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage[ref]{overcite}
\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}


\section*{Supplemental Materials}
Canonical ensemble (NVT) molecular dynamics calculations were
performed using the OOPSE molecular mechanics program.\cite{Meineke05}
All molecules were treated as rigid bodies, with orientational motion
propagated using the symplectic DLM integration method.  Details about
the implementation of this technique can be found in a recent
publication.\cite{Dullweber1997}

Thermodynamic integration is an established technique for
determination of free energies of condensed phases of
materials.\cite{Frenkel84,Hermens88,Meijer90,Baez95a,Vlot99}. This
method, implemented in the same manner by B\`{a}ez and Clancy, was
utilized to calculate the free energy of several ice crystals at 200 K
using the TIP3P, TIP4P, TIP5P, SPC/E, SSD/RF, and SSD/E water
models.\cite{Baez95a} Liquid state free energies at 300 and 400 K for
all of these water models were also determined using the same
technique in order to determine melting points and to generate phase
diagrams.  System sizes were 648 or 1728 molecules for ice B, 1024 or
1280 molecules for ice $I_h$, 1000 molecules for ice $I_c$, and 1024
molecules for both Ice-{\it i} and the liquid state simulations.  The
larger crystal sizes were necessary for simulations involving larger
cutoff values.  All simulations were carried out at densities which
correspond to a pressure of approximately 1 atm at their respective
temperatures.

Thermodynamic integration was utilized to calculate the Helmholtz free
energies ($A$) of the listed water models at various state points
using the OOPSE molecular dynamics program.\cite{Meineke05} This
method uses a sequence of simulations during which the system of
interest is converted into a reference system for which the free
energy is known analytically ($A_0$).  The difference in potential
energy between the reference system and the system of interest
($\Delta V$) is then integrated in order to determine the free energy
difference between the two states:
\begin{equation}
 A = A_0 + \int_0^1 \left\langle \Delta V \right\rangle_\lambda d\lambda.
\end{equation}
Simulations were {\it not} distributed uniformly along this path in
order to sufficiently sample the regions of greatest change in the
potential difference.  Typical integrations in this study consisted of
$\sim$25 simulations ranging from 300 ps (for the unaltered system) to
75 ps (near the reference state) in length.

For the thermodynamic integration of molecular crystals, the Einstein
crystal was chosen as the reference system.  In an Einstein crystal,
the molecules are restrained at their ideal lattice locations and
orientations. Using harmonic restraints, as applied by B\`{a}ez and
Clancy, the total potential for this reference crystal
($V_\mathrm{EC}$) is the sum of all the harmonic restraints,
\begin{equation}
V_\mathrm{EC} = \frac{K_\mathrm{v}r^2}{2} + \frac{K_\theta\theta^2}{2} +
\frac{K_\omega\omega^2}{2},
\end{equation}
where $K_\mathrm{v}$, $K_\mathrm{\theta}$, and $K_\mathrm{\omega}$ are
the spring constants restraining translational motion and deflection
of and rotation around the principle axis of the molecule
respectively.  These spring constants are typically calculated from
the mean-square displacements of water molecules in an unrestrained
ice crystal at 200 K.  For these studies, $K_\mathrm{r} = 4.29$ kcal
mol$^{-1}$, $K_\theta\ = 13.88$ kcal mol$^{-1}$, and $K_\omega\ =
17.75$ kcal mol$^{-1}$.  It is clear from Fig. \ref{waterSpring} that
the values of $\theta$ range from $0$ to $\pi$, while $\omega$ ranges
from $-\pi$ to $\pi$.  The partition function for a molecular crystal
restrained in this fashion can be evaluated analytically, and the
Helmholtz Free Energy ({\it A}) is given by
\begin{eqnarray}
A = E_m\ -\ kT\ln \left (\frac{kT}{h\nu}\right )^3&-&kT\ln \left
[\pi^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{A}kT}{h^2}\right
)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{B}kT}{h^2}\right
)^\frac{1}{2}\left (\frac{8\pi^2I_\mathrm{C}kT}{h^2}\right
)^\frac{1}{2}\right ] \nonumber \\ &-& kT\ln \left [\frac{kT}{2(\pi
K_\omega K_\theta)^{\frac{1}{2}}}\exp\left
(-\frac{kT}{2K_\theta}\right)\int_0^{\left (\frac{kT}{2K_\theta}\right
)^\frac{1}{2}}\exp(t^2)\mathrm{d}t\right ],
\label{ecFreeEnergy}
\end{eqnarray}
where $2\pi\nu = (K_\mathrm{v}/m)^{1/2}$, and $E_m$ is the minimum
potential energy of the ideal crystal.\cite{Baez95a}

\begin{figure}
\centering
\includegraphics[width=4in]{rotSpring.eps}
\caption{Possible orientational motions for a restrained molecule. 
$\theta$ angles correspond to displacement from the body-frame {\it
z}-axis, while $\omega$ angles correspond to rotation about the
body-frame {\it z}-axis.  $K_\theta$ and $K_\omega$ are spring
constants for the harmonic springs restraining motion in the $\theta$
and $\omega$ directions.}
\label{waterSpring}
\end{figure}

In the case of molecular liquids, the ideal vapor is chosen as the
target reference state.  There are several examples of liquid state
free energy calculations of water models present in the
literature.\cite{Hermens88,Quintana92,Mezei92,Baez95b} These methods
typically differ in regard to the path taken for switching off the
interaction potential to convert the system to an ideal gas of water
molecules.  In this study, we applied one of the most convenient
methods and integrated over the $\lambda^4$ path, where all
interaction parameters are scaled equally by this transformation
parameter.  This method has been shown to be reversible and provide
results in excellent agreement with other established
methods.\cite{Baez95b}

Near the cutoff radius ($0.85 * r_{cut}$), charge, dipole, and
Lennard-Jones interactions were gradually reduced by a cubic switching
function.  By applying this function, these interactions are smoothly
truncated, thereby avoiding the poor energy conservation which results
from harsher truncation schemes.  The effect of a long-range
correction was also investigated on select model systems in a variety
of manners.  For the SSD/RF model, a reaction field with a fixed
dielectric constant of 80 was applied in all
simulations.\cite{Onsager36} For a series of the least computationally
expensive models (SSD/E, SSD/RF, TIP3P, and SPC/E), simulations were
performed with longer cutoffs of 10.5, 12, 13.5, and 15 \AA\ to
compare with the 9 \AA\ cutoff results.  Finally, the effects of using
the Ewald summation were estimated for TIP3P and SPC/E by performing
single configuration Particle-Mesh Ewald (PME)
calculations~\cite{Tinker} for each of the ice polymorphs.  The
calculated energy difference in the presence and absence of PME was
applied to the previous results in order to predict changes to the
free energy landscape.

Additionally, $g_{OO}(r)$ and $S(\vec{q})$ plots were generated for
the two Ice-{\it i} variants (along with example ice $I_h$ and $I_c$
plots) at 77K, and they are shown in figures \ref{fig:gofr} and
\ref{fig:sofq}.  The $S(\vec{q})$ is related to a three dimensional
Fourier transform of the radial distribution function, which
simplifies to the following expression:

\begin{equation}
S(q) = 1 + 4\pi\rho\int_{0}^{\infty} r^2 \frac{\sin kr}{kr}g_{OO}(r)dr,
\label{sofqEq}
\end{equation}

where $\rho$ is the number density.  To obtain a good estimation of
$S(\vec{q})$, $g_{OO}(r)$ needs to extend to large $r$ values.  Thus,
simulations to obtain these plots were run for crystals eight times
the size of those used in the thermodynamic integrations.

\begin{figure}
\centering
\includegraphics[width=4in]{iceGofr.eps}
\caption{Radial distribution functions of ice $I_h$, $I_c$, and
Ice-{\it i} calculated from from simulations of the SSD/RF water model
at 77 K.  The Ice-{\it i} distribution function was obtained from
simulations composed of TIP4P water.}
\label{fig:gofr}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=4in]{sofq.eps}
\caption{Predicted structure factors for ice $I_h$, $I_c$, Ice-{\it i},
 and Ice-{\it i}$^\prime$ at 77 K.  The raw structure factors have
 been convoluted with a gaussian instrument function (0.075 \AA$^{-1}$
 width) to compensate for the trunction effects in our finite size
 simulations.}
\label{fig:sofq}
\end{figure}

\newpage

\bibliographystyle{jcp}
\bibliography{iceiSupplemental}

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