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}

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}
\includegraphics[width=\linewidth]{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 of 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 affects 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}
\includegraphics[width=\linewidth]{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}
\includegraphics[width=\linewidth]{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:	1879
Committed:	Thu Dec 9 23:15:07 2004 UTC (20 years, 9 months ago) by chrisfen
Content type:	application/x-tex
File size:	8648 byte(s)
Log Message:	improvements
#	User	Rev	Content
1	chrisfen	1872	%\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	chrisfen	1879	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	chrisfen	1872	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	chrisfen	1879	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	chrisfen	1872	\newpage
179
180			\bibliographystyle{jcp}
181			\bibliography{iceiSupplemental}
182
183
184			\end{document}