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}

\LARGE
Supplemental Material
\normalsize
\vspace{1cm}

Canonical ensemble (NVT) molecular dynamics calculations were
performed using the OOPSE molecular mechanics package.\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 illustrated 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
this 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 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 involves a sequence of simulations during
which the system of interest is converted into a reference system for
which the free energy is known analytically.  This transformation path
is then integrated in order to determine the free energy difference
between the two states:
\begin{equation}
\Delta A = \int_0^1\left\langle\frac{\partial V(\lambda 
)}{\partial\lambda}\right\rangle_\lambda d\lambda,
\end{equation}
where $V$ is the interaction potential and $\lambda$ is the
transformation parameter that scales the overall potential.
Simulations are distributed strategically along this path in order to
sufficiently sample the regions of greatest change in the potential.
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}
\begin{center}
\includegraphics[width=3in]{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{center}
\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}

Charge, dipole, and Lennard-Jones interactions were modified by a
cubic switching between 100\% and 85\% of the cutoff value (9 \AA).
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 provided by an Ewald summation were estimated for
TIP3P and SPC/E by performing single configuration calculations with
Particle-Mesh Ewald (PME) in the TINKER molecular mechanics software
package.\cite{Tinker} 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 on crystals eight times the
size of those used in the thermodynamic integrations.

\begin{figure}
\begin{center}
\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{center}
\end{figure}

\begin{figure}
\begin{center}
\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{center}
\end{figure}


\newpage

\bibliographystyle{jcp}
\bibliography{iceiSupplemental}


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