trunk/multipoleSFPaper/multipoleSFPaper.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
\documentclass[11pt]{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
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions
\let\Caption\caption
\renewcommand\caption[1]{%
        \Caption[#1]{}%
}
\setlength{\abovecaptionskip}{1.2in}
\setlength{\belowcaptionskip}{1.2in}

\begin{document}

\title{Pairwise Alternatives to the Ewald Sum: Applications
and Extension to Point Multipoles}

\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}
The damped, shifted-force electrostatic potential has been shown to
give nearly quantitative agreement with smooth particle mesh Ewald for
energy differences between configurations as well as for atomic force
and molecular torque vectors.  In this paper, we extend this technique
to handle interactions between electrostatic multipoles.  We also
investigate the effects of damped and shifted electrostatics on the
static, thermodynamic, and dynamic properties of liquid water and
various polymorphs of ice.  We provide a way of choosing the optimal
damping strength for a given cutoff radius that reproduces the static
dielectric constant in a variety of water models.
\end{abstract}

%\narrowtext 

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

\section{Introduction}

Over the past several years, there has been increasing interest in
pairwise methods for correcting electrostatic interactions in computer
simulations of condensed molecular
systems.\cite{Wolf99,Zahn02,Kast03,Beck05,Ma05,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} Wolf's method of cutoff neutralization is able to
obtain excellent agreement with Madelung energies in ionic
crystals.\cite{Wolf99} 

In a recent paper, we showed that simple modifications to Wolf's
method could give nearly quantitative agreement with smooth particle
mesh Ewald (SPME) for quantities of interest in Monte Carlo
(i.e. configurational energy differences) and Molecular Dynamics
(i.e. atomic force and molecular torque vectors).\cite{Fennell06} We
described the undamped and damped shifted potential (SP) and shifted
force (SF) techniques.  The potential for the damped form of the 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 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}
(In these potentials and in all electrostatic quantities that follow,
the standard $4 \pi \epsilon_{0}$ has been omitted for clarity.)

The damped SF method is an improvement over the 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 as well as the force effect of the charges
within $R_\textrm{c}$.

To complete the charge neutralization procedure, a self-neutralization
term is included in the potential. This term is constant (as long as
the charges and cutoff radius do not change), and exists outside the
normal pair-loop.  It can be thought of as a contribution from a
charge opposite in sign, but equal in magnitude, to the central
charge, but which has been spread out over the surface of the 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} = \frac{1}{2 R_\textrm{c}} \sum_{i=1}^N q_i^{2},
\label{eq:selfTerm}
\end{equation}
while for the damped case it can be written as
\begin{equation}
V_\textrm{self} = \left(\frac{\alpha}{\sqrt{\pi}} 
        + \frac{\textrm{erfc}(\alpha
R_\textrm{c})}{2R_\textrm{c}}\right) \sum_{i=1}^N q_i^{2}.
\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}

The SF, SP, and Wolf methods operate by neutralizing the total charge
contained within the cutoff sphere surrounding each particle.  This is
accomplished by creating image charges on the surface of the cutoff
sphere for each pair interaction computed within the sphere.  The
damping function applied to the potential is also an important method
for accelerating convergence.  In the case of systems involving
electrostatic distributions of higher order than point charges, the
question remains: How will the shifting and damping need to be
modified in order to accommodate point multipoles?

\section{Electrostatic Damping for Point
Multipoles}\label{sec:dampingMultipoles}

To apply the SF method for systems involving point multipoles, we
consider separately the two techniques (shifting and damping) which
contribute to the effectiveness of the DSF potential.

As noted above, shifting the potential and forces is employed to
neutralize the total charge contained within each cutoff sphere;
however, in a system composed purely of point multipoles, each cutoff
sphere is already neutral, so shifting becomes unnecessary.

In a mixed system of charges and multipoles, the undamped SF potential
needs only to shift the force terms between charges and smoothly
truncate the multipolar interactions with a switching function.  The
switching function is required for energy conservation, because a
discontinuity will exist in both the potential and forces at
$R_\textrm{c}$ in the absence of shifting terms. 

To damp the SF potential for point multipoles, we need to incorporate
the complimentary error function term into the standard forms of the
multipolar potentials.  We can determine the necessary damping
functions by replacing $1/r$ with $\mathrm{erfc}(\alpha r)/r$ in the
multipole expansion.  This procedure quickly becomes quite complex
with ``two-center'' systems, like the dipole-dipole potential, and is
typically approached using spherical harmonics.\cite{Hirschfelder67} A
simpler method for determining damped multipolar interaction
potentials arises when we adopt the tensor formalism described by
Stone.\cite{Stone02}

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

As one would do with the multipolar expansion, we can replace $r^{-1}$
with $\mathrm{erfc}(\alpha r)/r$ to obtain damped forms of the
electrostatic potentials.  Equation \ref{eq:tensorRank2} generates a
damped charge-dipole potential,
\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}

Equation \ref{eq:tensorRank3} generates a 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}
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 Smith and Aguado and Madden
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}$. 
$c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts
respectively. The forces for the damped dipole-dipole interaction, are
obtained from the next higher tensor, $T_{\alpha \beta \gamma}$,
\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}
Using the tensor formalism, we can dampen higher order multipolar
interactions using the same effective damping function that we use for
charge-charge interactions.  This allows us to include multipoles in
simulations involving damped electrostatic interactions.

As a practical consideration, we note that the evaluation of the error
function inside a pair loop can become quite costly.  In practice, we
pre-compute the $c_n(r)$ functions over a grid of $r$ values and use
cubic spline interpolation to obtain estimates of these functions when
necessary.  Using this procedure, the computational cost of damped
electrostatics is essentially the same as the undamped case.

\section{Applications of Damped Shifted-Froce Electrostatics}

Our earlier work on the SF method showed that it can give nearly
quantitive agreement with SPME-derived configurational energy
differences.  The force and torque vectors in identical configurations
are also nearly equivalent under the damped SF potential and
SPME.\cite{Fennell06} Although these measures bode well for the
performance of the SF method in both Monte Carlo and Molecular
Dynamics simulations, it would be helpful to have direct comparisons
of structural, thermodynamic, and dynamic quantities.  To address
this, we performed a detailed analysis of a group of simulations
involving water models (both point charge and multipolar) under a
number of different simulation conditions.

To provide the most difficult test for the damped SF method, we have
chosen a model that has been optimized for use with Ewald sum, and
have compared the simulated properties to those computed via Ewald.
It is well known that water models parametrized for use with the Ewald
sum give calculated properties that are in better agreement with
experiment.\cite{vanderSpoel98,Horn04,Rick04} For these reasons, we
chose the TIP5P-E water model for our comparisons involving point
charges.\cite{Rick04}  

The soft sticky dipole (SSD) family of water models is the perfect
test case for the point-multipolar extension of damped electrostatics.
SSD 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 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 the damped
electrostatic technique to see how the calculated properties change.

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 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.  In developing TIP5P-E, 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 compare these same properties calculated
from TIP5P-E using the Ewald sum with TIP5P-E using the damped SF
technique.  Our comparisons include the 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}$.

\subsection{The Density Maximum of TIP5P-E}\label{sec:t5peDensity}

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 SF method with
varying cutoff radii and damping coefficients. 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 SF methods show the average one-sigma
uncertainty of the density measurement, and this uncertainty is the
same for all the SF curves.}
\label{fig:t5peDensities}
\end{figure}
Figure \ref{fig:t5peDensities} shows the densities calculated for
TIP5P-E using differing electrostatic corrections overlaid with the
experimental values.\cite{CRC80} The densities when using the SF
technique are close to, but 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 SF sums showing slightly higher densities; however, it is
clear that the choice of cutoff radius plays a much more important
role in the resulting densities.

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 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.

\subsection{Liquid Structure of TIP5P-E}\label{sec:t5peLiqStructure}

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 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 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 do indicate that
$g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic
correction.

\subsection{Thermodynamic Properties of TIP5P-E}\label{sec:t5peThermo}

In addition to the density and structual features of the liquid, 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 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_{NPT} - \langle V\rangle^2_{NPT})}
                {k_BT\langle V\rangle_{NPT}},
\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_{NPT} 
                        - \langle V\rangle_{NPT}\langle H\rangle_{NPT})}
                        {k_BT^2\langle V\rangle_{NPT}}.
\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 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=5.8in]{./figures/t5peThermo.pdf}
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation 
and the 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}

For all of the properties computed, the trends with temperature
obtained when using the Ewald summation are reproduced with the SF
technique. One noticeable difference between the properties calculated
using the two methods are the lower $\Delta H_\textrm{vap}$ values
when using 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, 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 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 SF without damping. 

This points to the onset of a more frustrated or glassy behavior for
the undamped and weakly-damped SF simulations of TIP5P-E at
temperatures below 250~K than is seen from the Ewald sum.  Undamped SF
electrostatics has a stronger contribution from nearby charges.
Damping these local interactions or using a reciprocal-space method
makes the water less sensitive to ordering on a short length scale.
We can recover nearly quantitative agreement with the Ewald results by
increasing the damping constant.

The final thermodynamic property displayed in figure
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
between the Ewald and SF methods (and with experiment). It is known
that the dielectric constant is dependent upon and is quite sensitive
to the imposed boundary conditions.\cite{Neumann80,Neumann83} This is
readily apparent in the converged $\epsilon$ values accumulated for
the 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'' ($\kappa$) and real-space cutoff values also
have a significant effect on the calculated static dielectric constant
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.

\subsection{Dielectric Constants for TIP5P-E and SSD/RF}\label{sec:t5peDielectric} 

In the previous section, we observed that the choice of damping
coefficient plays a major role in the calculated dielectric constant
for the SF method.  Similar damping parameter behavior was observed in
the long-time correlated motions of the NaCl crystal.\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 quite sensitive to the
choice of damping parameter.  Since $\alpha$ is an adjustable
parameter, it would be best 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 for TIP5P-E and for a point-dipolar water model
(SSD/RF).  To calculate the static dielectric constant, we performed
5~ns $NPT$ calculations on systems of 512 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{Static dielectric constants for the TIP5P-E and SSD/RF water models at 298~K and 1~atm as a function of damping coefficient $\alpha$ and 
cutoff radius $R_\textrm{c}$. The color scale ranges from blue (10) to red (100).}
\vspace{6pt}
\begin{tabular}{ lccccccccc }
\toprule
\toprule
& \multicolumn{4}{c}{TIP5P-E} & & \multicolumn{4}{c}{SSD/RF} \\
& \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} & & \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\
\cmidrule(lr){2-5}  \cmidrule(lr){7-10}
$\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 & & 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.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.892222222222222, 0.5} & \cellcolor[rgb]{1, 0.704444444444444, 0.5} & \cellcolor[rgb]{1, 0.72, 0.5} & \cellcolor[rgb]{1, 0.6666666666667, 0.5}                                              & & \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.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       & & \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]{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       & & \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.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        & & \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, 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     & & \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.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.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.848888888888889, 1} & \cellcolor[rgb]{0.5, 0.973333333333333, 1} & \cellcolor[rgb]{0.5, 1, 0.793333333333333} & \cellcolor[rgb]{0.5, 1, 0.624444444444445}                               & & \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.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.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.696111111111111, 1} & \cellcolor[rgb]{0.5, 0.736333333333333, 1} & \cellcolor[rgb]{0.5, 0.775222222222222, 1} & \cellcolor[rgb]{0.5, 0.860666666666667, 1}                               & & \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.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   & & \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:DielectricMap}
\end{table}

The results of these calculations are displayed in table
\ref{tab:DielectricMap}.  The dielectric constants for both models
decrease linearly with increasing cutoff radii ($R_\textrm{c}$) and
increase linearly with increasing damping ($\alpha$).  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 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 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 SF technique, it might be interesting to
investigate the dielectric dependence of the real-space Ewald
parameters.

We have also mapped out the static dielectric constant of SSD/RF as a
function of $R_\textrm{c}$ and $\alpha$. 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.

Although it is tempting to choose damping parameters equivalent to the
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 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} 

We can use the same trend of $\alpha$ with $R_\textrm{c}$ for SSD/RF
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}  The static dielectric constants for SSD/RF
simulations with 9 and 12\AA\ $R_\textrm{c}$ values using their
respective damping parameters are 82.6$\pm$0.6 and 75$\pm$2.

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.

\subsection{Dynamic Properties of TIP5P-E}\label{sec:t5peDynamics}

To look at the dynamic properties of TIP5P-E when using the 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 from linear fits to the long-time
portion of the mean square displacement ($\langle r^{2}(t)
\rangle$).\cite{Allen87}  

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
axis $\alpha$.  The body-fixed reference frame used for our
orientational correlation functions has the $z$-axis running along the
HOH bisector, and the $y$-axis connecting the two hydrogen atoms.
$C_l^y$ is therefore calculated from the time evolution of a vector of
unit length pointing between the two hydrogen atoms.

From the orientation autocorrelation functions, we can obtain time
constants for rotational relaxation.  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=5.8in]{./figures/t5peDynamics.pdf}
\caption{Diffusion constants ({\it upper}) and reorientational time 
constants ({\it lower}) for TIP5P-E using the Ewald sum and 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 SF technique. The enhanced
diffusion (relative to experiment) at high temperatures are again the
product of the lower simulated densities and do not provide any
special insight into differences between the electrostatic summation
techniques. With the undamped SF technique, TIP5P-E tends to diffuse a
little faster than with the Ewald sum; however, the use of light to
moderate damping results in indistinguishable $D$ values. Though not
apparent in this figure, 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 phase with the SF technique at these lower temperatures, though
this change is more prominent with the {\it undamped} 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 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 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 SF results (25ps). A quick
calculation from a 10~ps trajectory with 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 lower 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 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{Comparison of Reaction Field and Damped Electrostatics for
SSD/RF} 

In section \ref{sec:dampingMultipoles}, we described a method for
applying the damped 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 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.


\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. FIXME] \\
& & $\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.

\subsection{Predictions of Ice Polymorph Stability}

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 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.  Since the
SF technique can be used as a simple and efficient replacement for the
Ewald summation, our final test of this method is to see if it is
capable of addressing the spurious stability of the Ice-$i$ phases
with the various common water models.  To this end, we have performed
thermodynamic integrations of all of the previously discussed ice
polymorphs using the 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 TIP5P and TIP4P 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 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.\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 different polymorphs to be nearly isoenergetic 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 using drastically different
charge distributions.  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 spurious 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
the solid forms of water than some of the other models.

\section{Conclusions}

This investigation of pairwise electrostatic summation techniques
shows that there is a viable and computationally efficient alternative
to the Ewald summation.  The SF method (eq. \ref{eq:DSFPot}) has
proven itself capable of reproducing structural, thermodynamic, and
dynamic quantities that are nearly quantitative matches to results
from far more expensive methods.  Additionally, we have now extended
the damping formalism to electrostatic multipoles, so the damped SF
potential can be used in systems that contain mixtures of charges and
point multipoles.

We have also provided a simple linear prescription for choosing
optimal damping parameters given a choice of cutoff radius.  The
damping parameters were chosen to obtain a static dielectric constant
as close as possible to the experimental value, which should be useful
for simulating the electrostatic screening properties of liquid water
accurately.  The linear formula for optimal damping was the same for
a complicated multipoint model as it was for a simple point-dipolar
model.  

As in all purely pairwise cutoff methods, the damped SF method is
expected to scale approximately {\it linearly} with system size, and
is easily parallelizable.  This should result in substantial
reductions in the computational cost of performing large simulations.
With the proper use of pre-computation and spline interpolation, the
damped SF method is exactly the same cost as a simple real-space
cutoff.

We are not suggesting that there is any flaw with the Ewald sum; in
fact, it is the standard by which the damped SF method has been
judged.  However, these results provide further evidence that in the
typical simulations performed today, the Ewald summation may no longer
be required to obtain the level of accuracy most researchers have come
to expect.

\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.  The authors would like
to thank Steve Corcelli and Ed Maginn for helpful discussions and
comments.

\newpage

\bibliographystyle{achemso}
\bibliography{multipoleSFPaper}


\end{document}
Revision:	3067
Committed:	Tue Oct 24 19:46:30 2006 UTC (18 years, 9 months ago) by gezelter
Content type:	application/x-tex
File size:	56149 byte(s)
Log Message:	more changes