trunk/multipoleSFPaper/multipoleSFPaper.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
\documentclass[12pt]{article}
\usepackage{tabls}
\usepackage{endfloat}
\usepackage[tbtags]{amsmath}
\usepackage{amsmath,bm}
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{colortbl}
\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}

\title{Pairwise electrostatic corrections for monopolar and multipolar systems}

\author{Christopher J. Fennell and J. Daniel Gezelter \\
Department of Chemistry and Biochemistry\\ University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle
%\doublespacing

\begin{abstract}
\end{abstract}

%\narrowtext 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

Over the past several years, there has been increased interest in
pairwise methods for correcting electrostatic interaction in computer
simulations of condensed molecular
systems.\cite{Wolf99,Zahn02,Kast03,Fennell06}

\section{Shifted Force Pairwise Electrostatic Correction}

In our previous look at pairwise electrostatic correction methods, we
described the undamped and damped shifted potential {\sc sp} and
shifted force {\sc sf} techniques.\cite{Fennell06} These techniques
were developed from the observations and efforts of Wolf {\it et al.}
and their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99}
The potential for the damped form of the {\sc sp} method is given by
\begin{equation}
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) 
\quad r\leqslant R_\textrm{c},
\label{eq:DSPPot}
\end{equation}
while the damped form of the {\sc sf} method is given by
\begin{equation}
\begin{split}
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
\frac{\mathrm{erfc}\left(\alpha r\right)}{r} 
- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ 
&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} 
+ \frac{2\alpha}{\pi^{1/2}}
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]} 
\quad r\leqslant R_\textrm{c}.
\label{eq:DSFPot}
\end{split}
\end{equation}
Overall, the damped {\sc sf} method is considered an improvement over
the {\sc sp} method because there is no discontinuity in the forces as
particles move out of the cutoff radius ($R_\textrm{c}$). This is
accomplished by shifting the forces (and potential) to zero at
$R_\textrm{c}$. This is analogous to neutralizing the charge and the
force effect of the charges within $R_\textrm{c}$. 

As we mentioned in the previous study, to complete the charge
neutralization procedure, a self-neutralization term needs to be
included in the potential. This term exists outside the pair-loop, and
can be thought of as a charge opposite in sign and equal in magnitude
to the central particle, spread out over the entire cutoff
sphere. This term is calculated via a single loop over all charges in
the system. For the undamped case, the self-term can be written as
\begin{equation}
V_\textrm{self} = \sum_{i=1}^N\frac{q_iq_i}{2R_\textrm{c}},
\label{eq:selfTerm}
\end{equation}
while for the damped case it can be written as
\begin{equation}
V_\textrm{Dself} = \sum_{i=1}^Nq_iq_i\left(\frac{\alpha}{\sqrt{\pi}} 
        + \frac{\textrm{erfc}(\alpha R_\textrm{c})}{2R_\textrm{c}}\right).
\label{eq:dampSelfTerm}
\end{equation}
The first term within the parentheses of equation
(\ref{eq:dampSelfTerm}) is identical to the self-term in the Ewald
summation, and comes from the utilization of the complimentary error
function for electrostatic damping.\cite{deLeeuw80,Wolf99}

\section{Multipolar Electrostatic Damping}\label{sec:dampingMultipoles}

As stated above, the {\sc sf} method operates by neutralizing the
charge pair force interactions (and potential) within the cutoff
sphere with shifting and by damping the electrostatic
interactions. Now we would like to consider an extension of these
techniques to include point multipole interactions. How will the
shifting and damping need to be modified in order to accommodate point
multipoles?

Of the two component techniques, the easiest to adapt is
shifting. Shifting is employed to neutralize the cutoff sphere;
however, in a system composed purely of point multipoles, the cutoff
sphere is already neutralized. This means that shifting is not
necessary between point multipoles. In a mixed system of monopoles and
multipoles, the undamped {\sc sf} potential needs only to shift the
force terms of the monopole and smoothly truncate the multipole
interactions with a switching function. The switching function is
required in order to conserve energy, because a discontinuity will
exist in both the potential and forces at $R_\textrm{c}$ in the
absence of shifting terms.

If we want to dampen the {\sc sf} potential, then we need to
incorporate the complimentary error function term into the multipole
potentials. The most direct way to do this is by replacing $r^{-1}$
with erfc$(\alpha r)\cdot r^{-1}$ in the multipole expansion. In the
multipole expansion, rather than considering only the interactions
between single point charges, the electrostatic interaction is
reformulated such that it describes the interaction between charge
distributions about central sites of the respective sets of
charges.\cite{Hirschfelder67} This procedure is what leads to the
familiar charge-dipole,
\begin{equation}
V_\textrm{cd} = -q_i\frac{\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}}{r^3_{ij}}
                = -q_i\mu_j\frac{\cos\theta}{r^2_{ij}}, 
\label{eq:chargeDipole}
\end{equation}
and dipole-dipole,
\begin{equation}
V_\textrm{dd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
                        (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}} -
                \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}},
\label{eq:dipoleDipole}
\end{equation}
interaction potentials.

Using the charge-dipole interaction as an example, if we insert
erfc$(\alpha r)\cdot r^{-1}$ in place of $r^{-1}$, a damped
charge-dipole results,
\begin{equation}
V_\textrm{Dcd} = -q_i\mu_j\frac{\cos\theta}{r^2_{ij}} c_1(r_{ij}),
\label{eq:dChargeDipole}
\end{equation}
where $c_1(r_{ij})$ is
\begin{equation}
c_1(r_{ij}) = \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}}
                + \textrm{erfc}(\alpha r_{ij}).
\label{eq:c1Func}
\end{equation}
Thus, $c_1(r_{ij})$ is the resulting damping term that modifies the
standard charge-dipole potential (Eq. (\ref{eq:chargeDipole})).  Note
that this damping term is dependent upon distance and not upon
orientation, and that it is acting on what was originally an
$r^{-3}$ function. By writing the damped form in this manner, we
can collect the damping into one function and apply it to the original
potential when damping is desired. This works well for potentials that
have only one $r^{-n}$ term (where $n$ is an odd positive integer);
but in the case of the dipole-dipole potential, there is one part
dependent on $r^{-3}$ and another dependent on $r^{-5}$. In order to
properly damping this potential, each of these parts is dampened with
separate damping functions. We can determine the necessary damping
functions by continuing with the multipole expansion; however, it
quickly becomes more complex with ``two-center'' systems, like the
dipole-dipole potential, and is typically approached with a spherical
harmonic formalism.\cite{Hirschfelder67} A simpler method for
determining these functions arises from adopting the tensor formalism
for expressing the electrostatic interactions.\cite{Stone02}

The tensor formalism for electrostatic interactions involves obtaining
the multipole interactions from successive gradients of the monopole
potential. Thus, tensors of rank one through three are
\begin{equation}
T = \frac{1}{4\pi\epsilon_0r_{ij}},
\label{eq:tensorRank1}
\end{equation}
\begin{equation}
T_\alpha = \frac{1}{4\pi\epsilon_0}\nabla_\alpha \frac{1}{r_{ij}},
\label{eq:tensorRank2}
\end{equation}
\begin{equation}
T_{\alpha\beta} = \frac{1}{4\pi\epsilon_0}
                        \nabla_\alpha\nabla_\beta \frac{1}{r_{ij}},
\label{eq:tensorRank3}
\end{equation}
where the form of the first tensor gives the monopole-monopole
potential, the second gives the monopole-dipole potential, and the
third gives the monopole-quadrupole and dipole-dipole
potentials.\cite{Stone02} Since the force is $-\nabla V$, the forces
for each potential come from the next higher tensor.

To obtain the damped electrostatic forms, we replace $r^{-1}$ with
erfc$(\alpha r)\cdot r^{-1}$. Equation \ref{eq:tensorRank2} generates
$c_1(r_{ij})$, just like the multipole expansion, while equation
\ref{eq:tensorRank3} generates $c_2(r_{ij})$, where
\begin{equation}
c_2(r_{ij}) = \frac{4\alpha^3r^3_{ij}e^{-\alpha^2r^2_{ij}}}{3\sqrt{\pi}}
                + \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}}
                + \textrm{erfc}(\alpha r_{ij}).
\end{equation} 
Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional
term. Continuing with higher rank tensors, we can obtain the damping
functions for higher multipole potentials and forces. Each subsequent
damping function includes one additional term, and we can simplify the
procedure for obtaining these terms by writing out the following
generating function,
\begin{equation}
c_n(r_{ij}) = \frac{2^n(\alpha r_{ij})^{2n-1}e^{-\alpha^2r^2_{ij}}}
                {(2n-1)!!\sqrt{\pi}} + c_{n-1}(r_{ij}),
\label{eq:dampingGeneratingFunc}
\end{equation}
where,
\begin{equation}
m!! \equiv \left\{ \begin{array}{l@{\quad\quad}l}
                m\cdot(m-2)\ldots 5\cdot3\cdot1 & m > 0\textrm{ odd} \\
                m\cdot(m-2)\ldots 6\cdot4\cdot2 & m > 0\textrm{ even} \\
                1 & m = -1\textrm{ or }0,
                \end{array}\right.
\label{eq:doubleFactorial}              
\end{equation}
and $c_0(r_{ij})$ is erfc$(\alpha r_{ij})$. This generating function
is similar in form to those obtained by researchers for the
application of the Ewald sum to
multipoles.\cite{Smith82,Smith98,Aguado03}

Returning to the dipole-dipole example, the potential consists of a
portion dependent upon $r^{-5}$ and another dependent upon
$r^{-3}$. In the damped dipole-dipole potential,
\begin{equation}
V_\textrm{Ddd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
                (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}}
                c_2(r_{ij}) -
                \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}}
                c_1(r_{ij}),
\label{eq:dampDipoleDipole}
\end{equation}
$c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts
respectively. The forces for the damped dipole-dipole interaction,
\begin{equation}
\begin{split}
F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
        (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})\mathbf{r}_{ij}}{r^7_{ij}} 
        c_3(r_{ij})\\ &- 
3\frac{\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}\cdot\hat{\boldsymbol{\mu}}_j +
        \boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}\cdot\hat{\boldsymbol{\mu}}_i +
        \boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}}
        {r^5_{ij}} c_2(r_{ij}),
\end{split}
\label{eq:dampDipoleDipoleForces}
\end{equation}
rely on higher order damping functions because we perform another
gradient operation. In this manner, we can dampen higher order
multipolar interactions along with the monopole interactions, allowing
us to include multipoles in simulations involving damped electrostatic
interactions.

\section{Applications of Pairwise Electrostatic Corrections: Water}

In our earlier work, we performed a statistical comparison of both
energetics and dynamics of the {\sc sf} and {\sc sp} methods, showing
that these pairwise correction techniques are nearly equivalent to the
Ewald summation in the simulation of general condensed
systems.\cite{Fennell06} It would be beneficial to have some specific
examples to better illustrate the similarities and differences of the
pairwise (specifically {\sc sf}) and Ewald techniques. To address
this, we decided to perform a detailed analysis involving water
simulations. Water is a good test in that there is a great deal of
detailed experimental information as well as an extensive library of
results from a variety of theoretical models. In choosing a model for
comparing these correction techniques, it is important to consider a
model that has been optimized for use with corrected electrostatics so
that the calculated properties are in better agreement with
experiment.\cite{vanderSpoel98,Horn04,Rick04} For this reason, we
chose the TIP5P-E water model.\cite{Rick04}

\subsection{TIP5P-E}

The TIP5P-E water model is a variant of Mahoney and Jorgensen's
five-point transferable intermolecular potential (TIP5P) model for
water.\cite{Mahoney00} TIP5P was developed to reproduce the density
maximum anomaly present in liquid water near 4$^\circ$C. As with many
previous point charge water models (such as ST2, TIP3P, TIP4P, SPC,
and SPC/E), TIP5P was parametrized using a simple cutoff with no
long-range electrostatic
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87}
Without this correction, the pressure term on the central particle
from the surroundings is missing. When this correction is included,
systems of these particles expand to compensate for this added
pressure term and under-predict the density of water under standard
conditions. When using any form of long-range electrostatic
correction, it has become common practice to develop or utilize a
reparametrized water model that corrects for this
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows
this practice and was optimized for use with the Ewald
summation.\cite{Rick04} In his publication, Rick preserved the
geometry and point charge magnitudes in TIP5P and focused on altering
the Lennard-Jones parameters to correct the density at 298~K. With the
density corrected, he compared common water properties for TIP5P-E
using the Ewald sum with TIP5P using a 9~\AA\ cutoff.

In the following sections, we compared these same water properties
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the
{\sc sf} technique.  In the above evaluation of the pairwise
techniques, we observed some flexibility in the choice of parameters.
Because of this, the following comparisons include the {\sc sf}
technique with a 12~\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and
0.2~\AA$^{-1}$, as well as a 9~\AA\ cutoff with an $\alpha$ =
0.2~\AA$^{-1}$.

\subsubsection{Density}\label{sec:t5peDensity}

As stated previously, the property that prompted the development of
TIP5P-E was the density at 1 atm.  The density depends upon the
internal pressure of the system in the $NPT$ ensemble, and the
calculation of the pressure includes a components from both the
kinetic energy and the virial. More specifically, the instantaneous
molecular pressure ($p(t)$) is given by
\begin{equation}
p(t) =  \frac{1}{\textrm{d}V}\sum_\mu
        \left[\frac{\mathbf{P}_{\mu}^2}{M_{\mu}} 
        + \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right],
\label{eq:MolecularPressure}
\end{equation}
where d is the dimensionality of the system, $V$ is the volume,
$\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$
is the position of the center of mass ($M_\mu$) of molecule $\mu$, and
$\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule
$\mu$.\cite{Melchionna93} The virial term (the right term in the
brackets of equation \ref{eq:MolecularPressure}) is directly dependent
on the interatomic forces.  Since the {\sc sp} method does not modify
the forces, the pressure using {\sc sp} will be identical to that
obtained without an electrostatic correction.  The {\sc sf} method
does alter the virial component and, by way of the modified pressures,
should provide densities more in line with those obtained using the
Ewald summation.

To compare densities, $NPT$ simulations were performed with the same
temperatures as those selected by Rick in his Ewald summation
simulations.\cite{Rick04} In order to improve statistics around the
density maximum, 3~ns trajectories were accumulated at 0, 12.5, and
25$^\circ$C, while 2~ns trajectories were obtained at all other
temperatures. The average densities were calculated from the later
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's
method for accumulating statistics, these sequences were spliced into
200 segments, each providing an average density. These 200 density
values were used to calculate the average and standard deviation of
the density at each temperature.\cite{Mahoney00}

\begin{figure}
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf}
\caption{Density versus temperature for the TIP5P-E water model when 
using the Ewald summation (Ref. \citen{Rick04}) and the {\sc sf} method
with various parameters. The pressure term from the image-charge shell
is larger than that provided by the reciprocal-space portion of the
Ewald summation, leading to slightly lower densities. This effect is
more visible with the 9~\AA\ cutoff, where the image charges exert a
greater force on the central particle. The error bars for the {\sc sf}
methods show the average one-sigma uncertainty of the density
measurement, and this uncertainty is the same for all the {\sc sf}
curves.}
\label{fig:t5peDensities}
\end{figure}
Figure \ref{fig:t5peDensities} shows the densities calculated for
TIP5P-E using differing electrostatic corrections overlaid on the
experimental values.\cite{CRC80} The densities when using the {\sc sf}
technique are close to, though typically lower than, those calculated
using the Ewald summation. These slightly reduced densities indicate
that the pressure component from the image charges at R$_\textrm{c}$
is larger than that exerted by the reciprocal-space portion of the
Ewald summation. Bringing the image charges closer to the central
particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than the
preferred 12~\AA\ R$_\textrm{c}$) increases the strength of the image
charge interactions on the central particle and results in a further
reduction of the densities.

Because the strength of the image charge interactions has a noticeable
effect on the density, we would expect the use of electrostatic
damping to also play a role in these calculations. Larger values of
$\alpha$ weaken the pair-interactions; and since electrostatic damping
is distance-dependent, force components from the image charges will be
reduced more than those from particles close the the central
charge. This effect is visible in figure \ref{fig:t5peDensities} with
the damped {\sc sf} sums showing slightly higher densities; however,
it is apparent that the choice of cutoff radius plays a much more
important role in the resulting densities.

As a final note, all of the above density calculations were performed
with systems of 512 water molecules. Rick observed a system size
dependence of the computed densities when using the Ewald summation,
most likely due to his tying of the convergence parameter to the box
dimensions.\cite{Rick04} For systems of 256 water molecules, the
calculated densities were on average 0.002 to 0.003~g/cm$^3$ lower. A
system size of 256 molecules would force the use of a shorter
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also
lower the densities. Moving to larger systems, as long as the
R$_\textrm{c}$ remains at a fixed value, we would expect the densities
to remain constant.

\subsubsection{Liquid Structure}\label{sec:t5peLiqStructure}

A common function considered when developing and comparing water
models is the oxygen-oxygen radial distribution function
($g_\textrm{OO}(r)$). The $g_\textrm{OO}(r)$ is the probability of
finding a pair of oxygen atoms some distance ($r$) apart relative to a
random distribution at the same density.\cite{Allen87} It is
calculated via
\begin{equation}
g_\textrm{OO}(r) = \frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i}
        \delta(\mathbf{r}-\mathbf{r}_{ij})\right\rangle,
\label{eq:GOOofR}
\end{equation}
where the double sum is over all $i$ and $j$ pairs of $N$ oxygen
atoms. The $g_\textrm{OO}(r)$ can be directly compared to X-ray or
neutron scattering experiments through the oxygen-oxygen structure
factor ($S_\textrm{OO}(k)$) by the following relationship:
\begin{equation}
S_\textrm{OO}(k) = 1 + 4\pi\rho
        \int_0^\infty r^2\frac{\sin kr}{kr}g_\textrm{OO}(r)\textrm{d}r.
\label{eq:SOOofK}
\end{equation}
Thus, $S_\textrm{OO}(k)$ is related to the Fourier-Laplace transform
of $g_\textrm{OO}(r)$.

The experimentally determined $g_\textrm{OO}(r)$ for liquid water has
been compared in great detail with the various common water models,
and TIP5P was found to be in better agreement than other rigid,
non-polarizable models.\cite{Sorenson00} This excellent agreement with
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To
check whether the choice of using the Ewald summation or the {\sc sf}
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at
298~K and 1~atm were determined for the systems compared in the
previous section.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf}
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298~K and 
1atm while using the Ewald summation (Ref. \cite{Rick04}) and the {\sc
sf} technique with varying parameters. Even with the reduced densities
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially
identical.}
\label{fig:t5peGofRs}
\end{figure}
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc
sf} technique with a various parameters are overlaid on the
$g_\textrm{OO}(r)$ while using the Ewald summation in
figure~\ref{fig:t5peGofRs}. The differences in density do not appear
to have any effect on the liquid structure as the $g_\textrm{OO}(r)$s
are indistinguishable. These results indicate that the
$g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic
correction.

\subsubsection{Thermodynamic Properties}\label{sec:t5peThermo}

In addition to the density, there are a variety of thermodynamic
quantities that can be calculated for water and compared directly to
experimental values. Some of these additional quantities include the
latent heat of vaporization ($\Delta H_\textrm{vap}$), the constant
pressure heat capacity ($C_p$), the isothermal compressibility
($\kappa_T$), the thermal expansivity ($\alpha_p$), and the static
dielectric constant ($\epsilon$). All of these properties were
calculated for TIP5P-E with the Ewald summation, so they provide a
good set for comparisons involving the {\sc sf} technique.

The $\Delta H_\textrm{vap}$ is the enthalpy change required to
transform one mole of substance from the liquid phase to the gas
phase.\cite{Berry00} In molecular simulations, this quantity can be
determined via
\begin{equation}
\begin{split}
\Delta H_\textrm{vap} &= H_\textrm{gas} - H_\textrm{liq.} \\
                      &= E_\textrm{gas} - E_\textrm{liq.} 
                        + p(V_\textrm{gas} - V_\textrm{liq.}) \\
                      &\approx -\frac{\langle U_\textrm{liq.}\rangle}{N}+ RT,
\end{split}
\label{eq:DeltaHVap}
\end{equation}
where $E$ is the total energy, $U$ is the potential energy, $p$ is the
pressure, $V$ is the volume, $N$ is the number of molecules, $R$ is
the gas constant, and $T$ is the temperature.\cite{Jorgensen98b} As
seen in the last line of equation (\ref{eq:DeltaHVap}), we can
approximate $\Delta H_\textrm{vap}$ by using the ideal gas for the gas
state. This allows us to cancel the kinetic energy terms, leaving only
the $U_\textrm{liq.}$ term. Additionally, the $pV$ term for the gas is
several orders of magnitude larger than that of the liquid, so we can
neglect the liquid $pV$ term.

The remaining thermodynamic properties can all be calculated from
fluctuations of the enthalpy, volume, and system dipole
moment.\cite{Allen87} $C_p$ can be calculated from fluctuations in the
enthalpy in constant pressure simulations via
\begin{equation}
\begin{split}
C_p     = \left(\frac{\partial H}{\partial T}\right)_{N,p} 
        = \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2},
\end{split}
\label{eq:Cp}
\end{equation}
where $k_B$ is Boltzmann's constant.  $\kappa_T$ can be calculated via
\begin{equation}
\begin{split}
\kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T}
         = \frac{(\langle V^2\rangle_{N,P,T} - \langle V\rangle^2_{N,P,T})}
                {k_BT\langle V\rangle_{N,P,T}},
\end{split}
\label{eq:kappa}
\end{equation}
and $\alpha_p$ can be calculated via
\begin{equation}
\begin{split}
\alpha_p        = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P}
                = \frac{(\langle VH\rangle_{N,P,T} 
                        - \langle V\rangle_{N,P,T}\langle H\rangle_{N,P,T})}
                        {k_BT^2\langle V\rangle_{N,P,T}}.
\end{split}
\label{eq:alpha}
\end{equation}
Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can
be calculated for systems of non-polarizable substances via
\begin{equation}
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV},
\label{eq:staticDielectric}
\end{equation}
where $\epsilon_0$ is the permittivity of free space and $\langle
M^2\rangle$ is the fluctuation of the system dipole
moment.\cite{Allen87} The numerator in the fractional term in equation
(\ref{eq:staticDielectric}) is the fluctuation of the simulation-box
dipole moment, identical to the quantity calculated in the
finite-system Kirkwood $g$ factor ($G_k$):
\begin{equation}
G_k = \frac{\langle M^2\rangle}{N\mu^2},
\label{eq:KirkwoodFactor}
\end{equation}
where $\mu$ is the dipole moment of a single molecule of the
homogeneous system.\cite{Neumann83,Neumann84,Neumann85} The
fluctuation term in both equation (\ref{eq:staticDielectric}) and
\ref{eq:KirkwoodFactor} is calculated as follows,
\begin{equation}
\begin{split}
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle 
                        - \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\
                   &= \langle M_x^2+M_y^2+M_z^2\rangle
                        - (\langle M_x\rangle^2 + \langle M_x\rangle^2 
                                + \langle M_x\rangle^2).        
\end{split}
\label{eq:fluctBoxDipole}
\end{equation}
This fluctuation term can be accumulated during the simulation;
however, it converges rather slowly, thus requiring multi-nanosecond
simulation times.\cite{Horn04} In the case of tin-foil boundary
conditions, the dielectric/surface term of the Ewald summation is
equal to zero. Since the {\sc sf} method also lacks this
dielectric/surface term, equation (\ref{eq:staticDielectric}) is still
valid for determining static dielectric constants.

All of the above properties were calculated from the same trajectories
used to determine the densities in section \ref{sec:t5peDensity}
except for the static dielectric constants. The $\epsilon$ values were
accumulated from 2~ns $NVE$ ensemble trajectories with system densities
fixed at the average values from the $NPT$ simulations at each of the
temperatures. The resulting values are displayed in figure
\ref{fig:t5peThermo}.
\begin{figure}
\centering
\includegraphics[width=4.5in]{./figures/t5peThermo.pdf}
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation 
and the {\sc sf} techniques along with the experimental values. Units
for the properties are kcal mol$^{-1}$ for $\Delta H_\textrm{vap}$,
cal mol$^{-1}$ K$^{-1}$ for $C_p$, 10$^6$ atm$^{-1}$ for $\kappa_T$,
and 10$^5$ K$^{-1}$ for $\alpha_p$. Ewald summation results are from
reference \cite{Rick04}. Experimental values for $\Delta
H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference
\cite{Kell75}. Experimental values for $C_p$ are from reference
\cite{Wagner02}. Experimental values for $\epsilon$ are from reference
\cite{Malmberg56}.}
\label{fig:t5peThermo}
\end{figure}

As observed for the density in section \ref{sec:t5peDensity}, the
property trends with temperature seen when using the Ewald summation
are reproduced with the {\sc sf} technique. One noticeable difference
between the properties calculated using the two methods are the lower
$\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be
expected due to the direct weakening of the electrostatic interaction
through forced neutralization. This results in an increase of the
intermolecular potential producing lower values from equation
(\ref{eq:DeltaHVap}). The slopes of these values with temperature are
similar to that seen using the Ewald summation; however, they are both
steeper than the experimental trend, indirectly resulting in the
inflated $C_p$ values at all temperatures.

Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ values
all overlap within error. As indicated for the $\Delta H_\textrm{vap}$
and $C_p$ results discussed in the previous paragraph, the deviations
between experiment and simulation in this region are not the fault of
the electrostatic summation methods but are due to the geometry and
parameters of the TIP5P class of water models. Like most rigid,
non-polarizable, point-charge water models, the density decreases with
temperature at a much faster rate than experiment (see figure
\ref{fig:t5peDensities}). This reduced density leads to the inflated
compressibility and expansivity values at higher temperatures seen
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability
and many-body effects are required in order for water models to
overcome differences between simulation-based and experimentally
determined densities at these higher
temperatures.\cite{Laasonen93,Donchev06}

At temperatures below the freezing point for experimental water, the
differences between {\sc sf} and the Ewald summation results are more
apparent. The larger $C_p$ and lower $\alpha_p$ values in this region
indicate a more pronounced transition in the supercooled regime,
particularly in the case of {\sc sf} without damping. This points to
the onset of a more frustrated or glassy behavior for TIP5P-E at
temperatures below 250~K in the {\sc sf} simulations, indicating that
disorder in the reciprocal-space term of the Ewald summation might act
to loosen up the local structure more than the image-charges in {\sc
sf}. The damped {\sc sf} actually makes a better comparison with
experiment in this region, particularly for the $\alpha_p$ values. The
local interactions in the undamped {\sc sf} technique appear to be too
strong since the property change is much more dramatic than the damped
forms, while the Ewald summation appears to weight the
reciprocal-space interactions at the expense the local interactions,
disagreeing with the experimental results. This observation is
explored further in section \ref{sec:t5peDynamics}.

The final thermodynamic property displayed in figure
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
between the Ewald summation and the {\sc sf} technique (and experiment
for that matter). It is known that the dielectric constant is
dependent upon and quite sensitive to the imposed boundary
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the
converged $\epsilon$ values accumulated for the {\sc sf}
simulations. Lack of a damping function results in dielectric
constants significantly smaller than those obtained using the Ewald
sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the
agreement considerably. It should be noted that the choice of the
``Ewald coefficient'' value also has a significant effect on the
calculated value when using the Ewald summation. In the simulations of
TIP5P-E with the Ewald sum, this screening parameter was tethered to
the simulation box size (as was the $R_\textrm{c}$).\cite{Rick04} In
general, systems with larger screening parameters reported larger
dielectric constant values, the same behavior we see here with {\sc
sf}; however, the choice of cutoff radius also plays an important
role. This connection is further explored below as optimal damping
coefficients for different choices of $R_\textrm{c}$ are determined
for {\sc sf} in order to best capture the dielectric behavior.

\subsubsection{Dielectric Constant}\label{sec:t5peDielectric}

In the previous, we observed that the choice of damping coefficient
plays a major role in the calculated dielectric constant.  This is not
too surprising given the results for damping parameter influence on
the long-time correlated motions of the NaCl crystal in our previous
study.\cite{Fennell06} The static dielectric constant is calculated
from the long-time fluctuations of the system's accumulated dipole
moment (Eq. (\ref{eq:staticDielectric})), so it is going to be quite
sensitive to the choice of damping parameter.  We would like to choose
optimal damping constants such that any arbitrary choice of cutoff
radius will properly capture the dielectric behavior of the liquid.

In order to find these optimal values, we mapped out the static
dielectric constant as a function of both the damping parameter and
cutoff radius.  To calculate the static dielectric constant, we
performed 5~ns $NPT$ calculations on systems of 512 TIP5P-E water
molecules and averaged over the converged region (typically the later
2.5~ns) of equation (\ref{eq:staticDielectric}).  The selected cutoff
radii ranged from 9, 10, 11, and 12~\AA , and the damping parameter
values ranged from 0.1 to 0.35~\AA$^{-1}$.

\begin{table}
\centering
\caption{TIP5P-E dielectric constant as a function of $\alpha$ and 
$R_\textrm{c}$. The color scale ranges from blue (10) to red (100).}
\vspace{6pt}
\begin{tabular}{ lcccc }
\toprule
\toprule
& \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\
\cmidrule(lr){2-5}
$\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 \\
\midrule
0.35 & \cellcolor[rgb]{1, 0.788888888888889, 0.5} 87.0 & \cellcolor[rgb]{1, 0.695555555555555, 0.5} 91.2 & \cellcolor[rgb]{1, 0.717777777777778, 0.5} 90.2 & \cellcolor[rgb]{1, 0.686666666666667, 0.5} 91.6 \\ 
& \cellcolor[rgb]{1, 0.892222222222222, 0.5} & \cellcolor[rgb]{1, 0.704444444444444, 0.5} & \cellcolor[rgb]{1, 0.72, 0.5} & \cellcolor[rgb]{1, 0.6666666666667, 0.5} \\ 
0.3 & \cellcolor[rgb]{1, 0.995555555555556, 0.5} 77.7 & \cellcolor[rgb]{1, 0.713333333333333, 0.5} 90.4 & \cellcolor[rgb]{1, 0.713333333333333, 0.5} 90.4 & \cellcolor[rgb]{1, 0.646666666666667, 0.5} 93.4 \\ 
0.275 & \cellcolor[rgb]{0.653333333333333, 1, 0.5} 61.9 & \cellcolor[rgb]{1, 0.933333333333333, 0.5} 80.5 & \cellcolor[rgb]{1, 0.811111111111111, 0.5} 86.0 & \cellcolor[rgb]{1, 0.766666666666667, 0.5} 88 \\ 
0.25 & \cellcolor[rgb]{0.537777777777778, 1, 0.5} 56.7 & \cellcolor[rgb]{0.795555555555555, 1, 0.5} 68.3 & \cellcolor[rgb]{1, 0.966666666666667, 0.5} 79 & \cellcolor[rgb]{1, 0.704444444444445, 0.5} 90.8 \\ 
0.225 & \cellcolor[rgb]{0.5, 1, 0.768888888888889} 42.9 & \cellcolor[rgb]{0.566666666666667, 1, 0.5} 58.0 & \cellcolor[rgb]{0.693333333333333, 1, 0.5} 63.7 & \cellcolor[rgb]{1, 0.937777777777778, 0.5} 80.3 \\ 
0.2 & \cellcolor[rgb]{0.5, 0.973333333333333, 1} 31.3 & \cellcolor[rgb]{0.5, 1, 0.842222222222222} 39.6 & \cellcolor[rgb]{0.54, 1, 0.5} 56.8 & \cellcolor[rgb]{0.735555555555555, 1, 0.5} 65.6 \\ 
& \cellcolor[rgb]{0.5, 0.848888888888889, 1} & \cellcolor[rgb]{0.5, 0.973333333333333, 1} & \cellcolor[rgb]{0.5, 1, 0.793333333333333} & \cellcolor[rgb]{0.5, 1, 0.624444444444445} \\ 
0.15 & \cellcolor[rgb]{0.5, 0.724444444444444, 1} 20.1 & \cellcolor[rgb]{0.5, 0.788888888888889, 1} 23.0 & \cellcolor[rgb]{0.5, 0.873333333333333, 1} 26.8 & \cellcolor[rgb]{0.5, 1, 0.984444444444444} 33.2 \\ 
& \cellcolor[rgb]{0.5, 0.696111111111111, 1} & \cellcolor[rgb]{0.5, 0.736333333333333, 1} & \cellcolor[rgb]{0.5, 0.775222222222222, 1} & \cellcolor[rgb]{0.5, 0.860666666666667, 1} \\ 
0.1 & \cellcolor[rgb]{0.5, 0.667777777777778, 1} 17.55 & \cellcolor[rgb]{0.5, 0.683777777777778, 1} 18.27 & \cellcolor[rgb]{0.5, 0.677111111111111, 1} 17.97 & \cellcolor[rgb]{0.5, 0.705777777777778, 1} 19.26 \\ 
\bottomrule
\end{tabular}
\label{tab:tip5peDielectricMap}
\end{table}
The results of these calculations are displayed in table
\ref{tab:tip5peDielectricMap}. Coloring of the table is based on the 
numerical values within the cells and was added to better facilitate
interpretation of the displayed data and highlight trends in the
dielectric constant behavior. This makes it easy to see that the
dielectric constant is linear with respect to $\alpha$ and
$R_\textrm{c}$ in the low to moderate damping regions, and the slope
is 0.025~\AA$^{-1}\cdot R_\textrm{c}^{-1}$. Another point to note is
that choosing $\alpha$ and $R_\textrm{c}$ identical to those used with
the Ewald summation results in the same calculated dielectric
constant. As an example, in the paper outlining the development of
TIP5P-E, the real-space cutoff and Ewald coefficient were tethered to
the system size, and for a 512 molecule system are approximately
12~\AA\ and 0.25~\AA$^{-1}$ respectively.\cite{Rick04} These
parameters resulted in a dielectric constant of 92$\pm$14, while with
{\sc sf} these parameters give a dielectric constant of
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and
0.35~\AA$^{-1}$, and these parameters resulted in a dielectric
constant equal to 63$\pm$1.\cite{Horn04} Calculations using {\sc sf}
with these parameters and this water model give a dielectric constant
of 61$\pm$1. Since the dielectric constant is dependent on $\alpha$
and $R_\textrm{c}$ with the {\sc sf} technique, it might be
interesting to investigate the dielectric dependence of the real-space
Ewald parameters.

Although it is tempting to choose damping parameters equivalent to
these Ewald examples, the results of our previous work indicate that
values this high are destructive to both the energetics and
dynamics.\cite{Fennell06} Ideally, $\alpha$ should not exceed
0.3~\AA$^{-1}$ for any of the cutoff values in this range. If the
optimal damping parameter is chosen to be midway between 0.275 and
0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$) at the 9~\AA\ cutoff, then the
linear trend with $R_\textrm{c}$ will always keep $\alpha$ below
0.3~\AA$^{-1}$ for the studied cutoff radii. This linear progression
would give values of 0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for
cutoff radii of 9, 10, 11, and 12~\AA. Setting this to be the default
behavior for the damped {\sc sf} technique will result in consistent
dielectric behavior for these and other condensed molecular systems,
regardless of the chosen cutoff radius. The static dielectric
constants for TIP5P-E simulations with 9 and 12\AA\ $R_\textrm{c}$
values using their respective damping parameters are 76$\pm$1 and
75$\pm$2. These values are lower than the values reported for TIP5P-E
with the Ewald sum; however, they are more in line with the values
reported by Mahoney {\it et al.} for TIP5P while using a reaction
field (RF) with an infinite RF dielectric constant
(81.5$\pm$1.6).\cite{Mahoney00} As a final note, aside from a slight
lowering of $\Delta H_\textrm{vap}$, using these $\alpha$ values does
not alter the other other thermodynamic properties.

\subsubsection{Dynamic Properties}\label{sec:t5peDynamics}

To look at the dynamic properties of TIP5P-E when using the {\sc sf}
method, 200~ps $NVE$ simulations were performed for each temperature
at the average density reported by the $NPT$ simulations using 9 and
12~\AA\ $R_\textrm{c}$ values using the ideal $\alpha$ values
determined above (0.2875 and 0.2125~\AA$^{-1}$). The self-diffusion
constants ($D$) were calculated using the mean square displacement
(MSD) form of the Einstein relation,
\begin{equation}
D = \lim_{t\rightarrow\infty} 
    \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
\label{eq:MSD}
\end{equation}
where $t$ is time, and $\mathbf{r}_i$ is the position of particle
$i$.\cite{Allen87} 

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/waterFrame.pdf}
\caption{Body-fixed coordinate frame for a water molecule. The 
respective molecular principle axes point in the direction of the
labeled frame axes.}
\label{fig:waterFrame}
\end{figure}
In addition to translational diffusion, orientational relaxation times
were calculated for comparisons with the Ewald simulations and with
experiments. These values were determined from the same 200~ps $NVE$
trajectories used for translational diffusion by calculating the
orientational time correlation function,
\begin{equation}
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t)
                \cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle,
\label{eq:OrientCorr}
\end{equation}
where $P_l$ is the Legendre polynomial of order $l$ and
$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along
principle axis $\alpha$. The principle axis frame for these water
molecules is shown in figure \ref{fig:waterFrame}. As an example,
$C_l^y$ is calculated from the time evolution of the unit vector
connecting the two hydrogen atoms.

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/exampleOrientCorr.pdf}
\caption{Example plots of the orientational autocorrelation functions 
for the first and second Legendre polynomials. These curves show the
time decay of the unit vector along the $y$ principle axis.}
\label{fig:OrientCorr}
\end{figure}
From the orientation autocorrelation functions, we can obtain time
constants for rotational relaxation. Figure \ref{fig:OrientCorr} shows
some example plots of orientational autocorrelation functions for the
first and second Legendre polynomials. The relatively short time
portions (between 1 and 3~ps for water) of these curves can be fit to
an exponential decay to obtain these constants, and they are directly
comparable to water orientational relaxation times from nuclear
magnetic resonance (NMR). The relaxation constant obtained from
$C_2^y(t)$ is of particular interest because it describes the
relaxation of the principle axis connecting the hydrogen atoms. Thus,
$C_2^y(t)$ can be compared to the intermolecular portion of the
dipole-dipole relaxation from a proton NMR signal and should provide
the best estimate of the NMR relaxation time constant.\cite{Impey82}

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/t5peDynamics.pdf}
\caption{Diffusion constants ({\it upper}) and reorientational time 
constants ({\it lower}) for TIP5P-E using the Ewald sum and {\sc sf}
technique compared with experiment. Data at temperatures less than
0$^\circ$C were not included in the $\tau_2^y$ plot to allow for
easier comparisons in the more relevant temperature regime.}
\label{fig:t5peDynamics}
\end{figure}
Results for the diffusion constants and orientational relaxation times
are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
the Ewald sum are reproduced with the {\sc sf} technique. The enhanced
diffusion at high temperatures are again the product of the lower
densities in comparison with experiment and do not provide any special
insight into differences between the electrostatic summation
techniques. With the undamped {\sc sf} technique, TIP5P-E tends to
diffuse a little faster than with the Ewald sum; however, use of light
to moderate damping results in indistinguishable $D$ values. Though
not apparent in this figure, {\sc sf} values at the lowest temperature
are approximately twice as slow as $D$ values obtained using the Ewald
sum. These values support the observation from section
\ref{sec:t5peThermo} that there appeared to be a change to a more
glassy-like phase with the {\sc sf} technique at these lower
temperatures, though this change seems to be more prominent with the
{\it undamped} {\sc sf} method, which has stronger local pairwise
electrostatic interactions.

The $\tau_2^y$ results in the lower frame of figure
\ref{fig:t5peDynamics} show a much greater difference between the {\sc
sf} results and the Ewald results. At all temperatures shown, TIP5P-E
relaxes faster than experiment with the Ewald sum while tracking
experiment fairly well when using the {\sc sf} technique, independent
of the choice of damping constant. Their are several possible reasons
for this deviation between techniques. The Ewald results were
calculated using shorter (10ps) trajectories than the {\sc sf} results
(25ps). A quick calculation from a 10~ps trajectory with {\sc sf} with
an $\alpha$ of 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in
$\tau_2^y$, placing the result more in line with that obtained using
the Ewald sum. This example supports this explanation; however,
recomputing the results to meet a poorer statistical standard is
counter-productive. Assuming the Ewald results are not entirely the
product of poor statistics, differences in techniques to integrate the
orientational motion could also play a role. {\sc shake} is the most
commonly used technique for approximating rigid-body orientational
motion,\cite{Ryckaert77} whereas in {\sc oopse}, we maintain and
integrate the entire rotation matrix using the {\sc dlm}
method.\cite{Meineke05} Since {\sc shake} is an iterative constraint
technique, if the convergence tolerances are raised for increased
performance, error will accumulate in the orientational
motion. Finally, the Ewald results were calculated using the $NVT$
ensemble, while the $NVE$ ensemble was used for {\sc sf}
calculations. The additional mode of motion due to the thermostat will
alter the dynamics, resulting in differences between $NVT$ and $NVE$
results. These differences are increasingly noticeable as the
thermostat time constant decreases.


\subsection{SSD/RF}

In section \ref{sec:dampingMultipoles}, we described a method for
applying the damped {\sc sf} technique to systems containing point
multipoles. The soft sticky dipole (SSD) family of water models is the
perfect test case because of the dipole-dipole (and
charge-dipole/quadrupole) interactions that are present. As alluded to
in the name, soft sticky dipole water models are single point
molecules that consist of a ``soft'' Lennard-Jones sphere, a
point-dipole, and a tetrahedral function for capturing the hydrogen
bonding nature of water - a spherical harmonic term for water-water
tetrahedral interactions and a point-quadrupole for interactions with
surrounding charges. Detailed descriptions of these models can be
found in other studies.\cite{Liu96b,Chandra99,Tan03,Fennell04}

In deciding which version of the SSD model to use, we need only
consider that the {\sc sf} technique was presented as a pairwise
replacement for the Ewald summation. It has been suggested that models
parametrized for the Ewald summation (like TIP5P-E) would be
appropriate for use with a reaction field and vice versa.\cite{Rick04}
Therefore, we decided to test the SSD/RF water model, which was
parametrized for use with a reaction field, with this damped
electrostatic technique to see how the calculated properties change.

\subsubsection{Dipolar Damping}

\begin{table}
\caption{Properties of SSD/RF when using different electrostatic 
correction methods.}
\footnotesize
\centering
\begin{tabular}{ llccc } 
\toprule
\toprule
& & Reaction Field [Ref. \citen{Fennell04}] & Damped Electrostatics & Experiment [Ref.] \\
& & $\epsilon = 80$ & $R_\textrm{c} = 12$\AA ; $\alpha = 0.2125$~\AA$^{-1}$ & \\
\midrule
$\rho$ & (g cm$^{-3}$) & 0.997(1) & 1.004(4) & 0.997 [\citen{CRC80}]\\
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 23.8(2) & 27(1) & 18.005 [\citen{Wagner02}] \\
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.32(6) & 2.33(2) & 2.299 [\citen{Mills73}]\\
$n_C$ & & 4.4 & 4.2 & 4.7 [\citen{Hura00}]\\ 
$n_H$ & & 3.7 & 3.7 & 3.5 [\citen{Soper86}]\\
$\tau_1$ & (ps) & 7.2(4) & 5.86(8) & 5.7 [\citen{Eisenberg69}]\\
$\tau_2$ & (ps) & 3.2(2) & 2.45(7) & 2.3 [\citen{Krynicki66}]\\
\bottomrule
\end{tabular}
\label{tab:dampedSSDRF}
\end{table}
The properties shown in table \ref{tab:dampedSSDRF} compare
quite well. The average density shows a modest increase when
using damped electrostatics in place of the reaction field. This comes
about because we neglect the pressure effect due to the surroundings
outside of the cutoff, instead relying on screening effects to
neutralize electrostatic interactions at long distances. The $C_p$
also shows a slight increase, indicating greater fluctuation in the
enthalpy at constant pressure. The only other differences between the
damped and reaction field results are the dipole reorientational time
constants, $\tau_1$ and $\tau_2$. When using damped electrostatics,
the water molecules relax more quickly and exhibit behavior very
similar to the experimental values. These results indicate that not
only is it reasonable to use damped electrostatics with SSD/RF, it is
recommended if capturing realistic dynamics is of primary
importance. This is an encouraging result because the damping methods
are more generally applicable than reaction field. Using damping,
SSD/RF can be used effectively with mixed systems, such as dissolved
ions, dissolved organic molecules, or even proteins.

\subsubsection{Dielectric Constant}

\begin{table}
\centering
\caption{SSD/RF dielectric constant as a function of $\alpha$ and $R_\textrm{c}$. The color scale ranges from blue (10) to red (100).}
\vspace{6pt}
\begin{tabular}{ lcccc }
\toprule
\toprule
& \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\
\cmidrule(lr){2-5}
$\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 \\
\midrule
0.35 & \cellcolor[rgb]{1, 0.5, 0.5} 116 & \cellcolor[rgb]{1, 0.5, 0.5} 119.2 & \cellcolor[rgb]{1, 0.5, 0.5} 131.4 & \cellcolor[rgb]{1, 0.5, 0.5} 130 \\ 
& \cellcolor[rgb]{1, 0.5, 0.5} & \cellcolor[rgb]{1, 0.5, 0.5} & \cellcolor[rgb]{1, 0.5, 0.5} & \cellcolor[rgb]{1, 0.5, 0.5} \\ 
0.3 & \cellcolor[rgb]{1, 0.5, 0.5} 100 & \cellcolor[rgb]{1, 0.5, 0.5} 118.8 & \cellcolor[rgb]{1, 0.5, 0.5} 116 & \cellcolor[rgb]{1, 0.5, 0.5} 122 \\ 
0.275 & \cellcolor[rgb]{1, 1, 0.5} 77.5 & \cellcolor[rgb]{1, 0.5, 0.5} 105 & \cellcolor[rgb]{1, 0.5, 0.5} 118 & \cellcolor[rgb]{1, 0.5, 0.5} 125.2 \\ 
0.25 & \cellcolor[rgb]{0.5, 1, 0.582222222222222} 51.3 & \cellcolor[rgb]{1, 0.993333333333333, 0.5} 77.8 & \cellcolor[rgb]{1, 0.522222222222222, 0.5} 99 & \cellcolor[rgb]{1, 0.5, 0.5} 113 \\ 
0.225 & \cellcolor[rgb]{0.5, 0.984444444444444, 1} 31.8 & \cellcolor[rgb]{0.5, 1, 0.586666666666667} 51.1 & \cellcolor[rgb]{1, 0.995555555555556, 0.5} 77.7 & \cellcolor[rgb]{1, 0.5, 0.5} 108.1 \\ 
0.2 & \cellcolor[rgb]{0.5, 0.698666666666667, 1} 18.94 & \cellcolor[rgb]{0.5, 0.946666666666667, 1} 30.1 & \cellcolor[rgb]{0.5, 1, 0.704444444444445} 45.8 & \cellcolor[rgb]{0.893333333333333, 1, 0.5} 72.7 \\ 
& \cellcolor[rgb]{0.5, 0.599333333333333, 1} & \cellcolor[rgb]{0.5, 0.732666666666667, 1} & \cellcolor[rgb]{0.5, 0.942111111111111, 1} & \cellcolor[rgb]{0.5, 1, 0.695555555555556} \\ 
0.15 & \cellcolor[rgb]{0.5, 0.5, 1} 8.29 & \cellcolor[rgb]{0.5, 0.518666666666667, 1} 10.84 & \cellcolor[rgb]{0.5, 0.588666666666667, 1} 13.99 & \cellcolor[rgb]{0.5, 0.715555555555556, 1} 19.7 \\ 
& \cellcolor[rgb]{0.5, 0.5, 1} & \cellcolor[rgb]{0.5, 0.509333333333333, 1} & \cellcolor[rgb]{0.5, 0.544333333333333, 1} & \cellcolor[rgb]{0.5, 0.607777777777778, 1} \\ 
0.1 & \cellcolor[rgb]{0.5, 0.5, 1} 4.96 & \cellcolor[rgb]{0.5, 0.5, 1} 5.46 & \cellcolor[rgb]{0.5, 0.5, 1} 6.04 & \cellcolor[rgb]{0.5, 0.5, 1} 6.82 \\ 
\bottomrule
\end{tabular}
\label{tab:ssdrfDielectricMap}
\end{table}
In addition to the properties tabulated in table
\ref{tab:dampedSSDRF}, we mapped out the static dielectric constant of
SSD/RF as a function of $R_\textrm{c}$ and $\alpha$ in the same
fashion as in section \ref{sec:t5peDielectric} with TIP5P-E. It is
apparent from this table that electrostatic damping has a more
pronounced effect on the dipolar interactions of SSD/RF than the
monopolar interactions of TIP5P-E. The dielectric constant covers a
much wider range and has a steeper slope with increasing damping
parameter. We can use the same trend of $\alpha$ with $R_\textrm{c}$
used with TIP5P-E, and for a 12~\AA\ $R_\textrm{c}$, the resulting
dielectric constant is 82.6$\pm$0.6. This value compares very
favorably with the experimental value of 78.3.\cite{Malmberg56} This
is not surprising given that early studies of the SSD model indicate a
static dielectric constant around 81.\cite{Liu96}

\section{Application of Pairwise Electrostatic Corrections: Imaginary Ice}

In an earlier work, we performed a series of free energy calculations
on several ice polymorphs which are stable or metastable at low
pressures, one of which (Ice-$i$) we observed in spontaneous
crystallizations of an SSD type single point water
model.\cite{Fennell05} In this study, a distinct dependence of the
free energies on the interaction cutoff and correction technique was
observed. Being that the {\sc sf} technique can be used as a simple
and efficient replacement for the Ewald summation, it would be
interesting to apply it in addressing the question of the free
energies of these ice polymorphs with varying water models. To this
end, we set up thermodynamic integrations of all of the previously
discussed ice polymorphs using the {\sc sf} technique with a cutoff
radius of 12~\AA\ and an $\alpha$ of 0.2125~\AA . These calculations
were performed on TIP5P-E and TIP4P-Ew (variants of the root models
optimized for the Ewald summation) as well as SPC/E, and SSD/RF.

\begin{table}
\centering
\caption{Helmholtz free energies of ice polymorphs at 1~atm and 200~K 
using the damped {\sc sf} electrostatic correction method with a
variety of water models.}
\begin{tabular}{ lccccc }
\toprule
\toprule
Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\
\cmidrule(lr){2-6}
& \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\
\midrule 
TIP5P-E & -11.98(4) & -11.96(4) & -11.87(3) & - & -11.95(3) \\
TIP4P-Ew & -13.11(3) & -13.09(3) & -12.97(3) & - & -12.98(3) \\
SPC/E & -12.99(3) & -13.00(3) & -13.03(3) & - & -12.99(3) \\
SSD/RF & -11.83(3) & -11.66(4) & -12.32(3) & -12.39(3) & - \\
\bottomrule
\end{tabular}
\label{tab:dampedFreeEnergy}
\end{table}
The results of these calculations in table \ref{tab:dampedFreeEnergy}
show similar behavior to the Ewald results in the previous study, at
least for SSD/RF and SPC/E which are present in both.\cite{Fennell05}
The Helmholtz free energies of the ice polymorphs with SSD/RF order in
the same fashion, with Ice-$i$ having the lowest free energy; however,
the Ice-$i$ and ice B free energies are quite a bit closer (nearly
isoenergetic). The SPC/E results show the near isoenergetic behavior
when using the Ewald summation.\cite{Fennell05} Ice B has the lowest
Helmholtz free energy; however, all the polymorph results overlap
within error.

The most interesting results from these calculations come from the
more expensive TIP4P-Ew and TIP5P-E results. Both of these models were
optimized for use with an electrostatic correction and are
geometrically arranged to mimic water following two different
ideas. In TIP5P-E, the primary location for the negative charge in the
molecule is assigned to the lone-pairs of the oxygen, while TIP4P-Ew
places the negative charge near the center-of-mass along the H-O-H
bisector. There is some debate as to which is the proper choice for
the negative charge location, and this has in part led to a six-site
water model that balances both of these options.\cite{Vega05,Nada03}
The limited results in table \ref{tab:dampedFreeEnergy} support the
results of Vega {\it et al.}, which indicate the TIP4P charge location
geometry is more physically valid.\cite{Vega05} With the TIP4P-Ew
water model, the experimentally observed polymorph (ice
I$_\textrm{h}$) is the preferred form with ice I$_\textrm{c}$ slightly
higher in energy, though overlapping within error. Additionally, the
less realistic ice B and Ice-$i^\prime$ structures are destabilized
relative to these polymorphs. TIP5P-E shows similar behavior to SPC/E,
where there is no real free energy distinction between the various
polymorphs, because many overlap within error. While ice B is close in
free energy to the other polymorphs, these results fail to support the
findings of other researchers indicating the preferred form of TIP5P
at 1~atm is a structure similar to ice
B.\cite{Yamada02,Vega05,Abascal05} It should be noted that we are
looking at TIP5P-E rather than TIP5P, and the differences in the
Lennard-Jones parameters could be a reason for this dissimilarity.
Overall, these results indicate that TIP4P-Ew is a better mimic of
real water than these other models when studying crystallization and
solid forms of water.

\section{Conclusions}

\section{Acknowledgments}
Support for this project was provided by the National Science
Foundation under grant CHE-0134881. Computation time was provided by
the Notre Dame Center for Research Computing.

\newpage

\bibliographystyle{achemso}
\bibliography{multipoleSFPaper}


\end{document}
Revision:	3065
Committed:	Tue Oct 24 06:27:31 2006 UTC (19 years, 7 months ago) by chrisfen
Content type:	application/x-tex
File size:	56789 byte(s)
Log Message:	fixed the bib file and some references