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\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
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.  Additionally, 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}

\newpage

%\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, where $\alpha$ is the adjustable damping parameter, is given
by
\begin{equation}
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) 
\quad r_{ij}\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_{ij}\right)}{r_{ij}} 
- \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_{ij}-R_\mathrm{c}\right)\ \Biggr{]} 
\quad r_{ij}\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 added to 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, 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 shifting the potential functions to generate image
charges on the surface of the cutoff sphere for each pair interaction
computed within $R_\textrm{c}$.  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.  It would
still be possible, however, to subtract out the effects of image
multipoles.  This is essentially what is done for dipoles in the
reaction field methods, and the technique could be extended to higher
order multipoles.   

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 dampen 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_{ij}$ with $\mathrm{erfc}(\alpha r_{ij})/r_{ij}$ 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_{ij}^{-1}$
with $\mathrm{erfc}(\alpha r_{ij})/r_{ij}$ 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\frac{\mathbf{r}_{ij}\cdot\boldsymbol{\mu}_j}{r^3_{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
recurrence relation,
\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,\cite{Smith82,Smith98}
Toukmaji {\it et al.}\cite{Toukmaji00}, Aguado and
Madden,\cite{Aguado03},  and Sagui {\it et al.}\cite{Sagui04} for the 
application of the Ewald sum to multipoles.

Returning to the dipole-dipole example, the potential consists of a
portion dependent upon $r_{ij}^{-5}$ and another dependent upon
$r_{ij}^{-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\boldsymbol{\mu}_j +
        (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})\cdot\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.  In general,
if the multipolar potentials are left in $\mathbf{r}_{ij}/r_{ij}$
form, instead of reducing them to the more common angular forms which
use $\hat{\mathbf{r}}_{ij}$ (or the resultant angles), one may simply replace
any $1/r_{ij}^{2n+1}$ dependence with $c_n(r_{ij}) / r_{ij}^{2n+1}$ to
obtain the damped version of that multipolar potential.

As a practical consideration, we note that the evaluation of the
complementary 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 only marginally
more costly than the undamped case.

\section{Applications of Damped Shifted-Force 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 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,
the system expands to compensate for this added pressure and therefore
under-predicts 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}$.

Moving beyond point-charge electrostatics, 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.

\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 latter
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 representative error bar for the SF methods shows 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 than the undamped
case; 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 oxygen-oxygen pair correlation function
($g_\textrm{OO}(r)$) for liquid water has been compared in great
detail with predictions of 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, we calculated this correlation function at 298~K and
1~atm for the parameters used in the previous section.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf}
\caption{The oxygen-oxygen pair correlation functions calculated for
TIP5P-E at 298~K and 1~atm while using the Ewald summation
(Ref. \citen{Rick04}) and the SF technique with varying
parameters. Even with the lower densities obtained using the SF
technique, the correlation functions are essentially identical.}
\label{fig:t5peGofRs}
\end{figure}
The pair correlation functions calculated for TIP5P-E while using the
SF technique with various parameters are overlaid on the same function
obtained 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 correlation
functions 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 of reference data 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 \citen{Rick04}. Experimental values for $\Delta
H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference
\citen{Kell75}. Experimental values for $C_p$ are from reference
\citen{Wagner02}. Experimental values for $\epsilon$ are from reference
\citen{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 at these temperatures.  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 shorter 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 general, aqueous systems with larger screening
parameters will report larger values for the dielectric constant, the
same behavior we see here with SF.  Although the received wisdom
(CITATION) on the Ewald parameter is that large enough Fourier-space
grids can make up for the effects of overdamping the real-space
portion of the Ewald sum, this wisdom does not appear to be borne out
by actual numerical results.  To illustrate this point, we have
computed the static dielectric constants from 216-molecule SPC water
simulations (9 \AA\ cutoff radius, 3 ns data collection times, PME
with fourth-order spline to smooth the reciprocal space portion of the
sum).  We carried out these simulations with a range of different
Ewald coefficients.  For each value of the Ewald coefficient, we also
varied the number of Fourier vectors used to compute the $k$-space
sum.  At the largest Fourier-grid (48x48x48 Fourier-vectors), this
corresponds to a 0.39 \AA\ grid spacing.  The results for these
calculations are shown in figure \ref{fig:ewaldDielectric}.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/ewaldDielectric.pdf}
\caption{The dielectric constant of the SPC water model using the Ewald
summation as a function of ({\it left}) the Ewald coefficient
($\kappa$) and ({\it right}) the $k$-Space grid detail.  All
calculations were performed with a 216 water molecule system using a
fixed cutoff radius of 9 \AA\ at 300 K.}
\label{fig:ewaldDielectric}
\end{figure}

It is clear that the choice of the Ewald coefficient actually does
have a dramatic effect on the static dielectric of water, and that
increasing the Fourier grid to 3 times the typical value does not
bring these dielectric constants back into agreement with each other.
It should be noted that this behavior would be troubling in {\it any}
method for calculating electrostatic interactions.  Ideally, a method
should give fixed values for important properties like the static
dielectric and should be insensitive to the choice of convergence
parameters like the Ewald coefficient and $\alpha$.  From what we have
observed, neither the Ewald summation nor the real-space shifted
methods appear to be ``black boxes'' at this point.  Both require
careful consideration of the choice of parameters on the property
being studied.

\subsection{Optimal Damping Coefficients for Damped
Electrostatics}\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 latter 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 with increasing cutoff radii ($R_\textrm{c}$) and increase
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 dependence of the static dielectric
constant on the choice of convergence parameters ($R_\textrm{c}$ and
$\kappa$) utilized in most implementations of the Ewald sum.

It is also 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 to obtain quantitative agreement, 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. 

It is possible to use a simple energetic criterion for choosing an
optimal value of $\alpha$ given a cutoff radius ($R_\textrm{c}$) and
an energetic tolerance.  This is identical to how the Ewald
coefficient ($\kappa$) is chosen in many of the common simulation
packages (GROMACS,\cite{Gromacs} AMBER,\cite{AMBER} and
TINKER~\cite{TINKER}).  In the case of the damped method described in
this work, the value of the tolerance,
\begin{equation}
\textrm{tolerance} = \frac{V_\textrm{damped}(R_\textrm{c})}{V_\textrm{coulomb}(R_\textrm{c})}
= \textrm{erfc}(\alpha R_\textrm{c}),
\end{equation}
which does best at reproducing the static dielectric constant in a
variety of water models is $3.1~\times~10^{-4}$.  Using this
tolerance, good choices of $\alpha$ for cutoff radii of 9 and 12
\AA\ are 0.2834 and 0.2125 \AA$^{-1}$, respectively.  

As a final note on optimal damping parameters, aside from a slight
lowering of $\Delta H_\textrm{vap}$, using optimal $\alpha$ values
does not alter any of the 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 obtained from the $NPT$ simulations.  $R_\textrm{c}$
values of 9 and 12~\AA\ and the ideal $\alpha$ values determined above
(0.2834 and 0.2125~\AA$^{-1}$) were used for the damped
electrostatics. 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^\gamma(t)
                \cdot\hat{\mathbf{u}}_i^\gamma(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 $\gamma$.  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.  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 the SF simulations result in a
slightly more viscous supercooled region than is obtained using the
Ewald sum.

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). There are several possible reasons for
this deviation between techniques. The Ewald results were calculated
using shorter trajectories (10~ps) than the SF results (200~ps).
Calculation of these SF values from 10~ps trajectories (with
subsequently lower accuracy) showed a 0.4~ps drop in $\tau_2^y$,
placing the result more in line with that obtained using the Ewald
sum.  Recomputing correlation times to meet a lower statistical
standard is counter-productive, however.  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 motion due to the extended variable (the thermostat)
will always alter the dynamics, resulting in differences between $NVT$
and $NVE$ results. These differences are increasingly noticeable as
the time constant for the thermostat decreases.

\subsection{Comparison of Reaction Field and Damped Electrostatics for
SSD/RF} 

SSD/RF was parametrized for use with a reaction field, which is a
common and relatively inexpensive way of handling long-range
electrostatic corrections in dipolar systems.\cite{Onsager36}
Although there is no reason to expect that damped electrostatics will
behave in a similar fashion to the reaction field, it is well known
that model that are parametrized for use with Ewald do better than
unoptimized models under the influence of a reaction
field.\cite{Rick04}  We compared a number of properties of SSD/RF that
had previously been computed using a reaction field with those same
values under damped electrostatics. 

The properties shown in table \ref{tab:dampedSSDRF} show that using
damped electrostatics can result in even better agreement with
experiment than is obtained via reaction field.  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 electrostatics and the 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 closer to the experimental
values. These results indicate that since there is no need to specify
an external dielectric constant with the damped electrostatics, it is
almost certainly a better choice for dipolar simulations than the
reaction field method.  Additionally, by using damped electrostatics
instead of reaction field, SSD/RF can be used effectively with mixed
charge / dipolar systems, such as dissolved ions, dissolved organic
molecules, or even proteins.

\begin{table}
\caption{Properties of SSD/RF when using reaction field or damped
electrostatics.  Simulations were carried out at 298~K, 1~atm, and
with 512 molecules.}
\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}

\subsection{Predictions of Ice Polymorph Stability}

In an earlier paper, 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 early version of the SSD/RF water
model.\cite{Fennell05} In this study, a distinct dependence of the
computed free energies on the cutoff radius and electrostatic
summation method 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.  This is the same
behavior observed using an Ewald correction.\cite{Fennell05} Ice B has
the lowest Helmholtz free energy for SPC/E; however, all the polymorph
results overlap within the error estimates.

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 performs better
under some circumstances.\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. 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 were
looking at TIP5P-E rather than TIP5P, and the differences in the
Lennard-Jones parameters could cause this discrepancy.  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 (equation (\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 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 formula for optimal damping was the same for a
complicated multipoint model as it was for a simple point-dipolar
model, and is identical to energetic tolerance methods commonly used
to choose the Ewald coefficient. 

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 essentially 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:	3085
Committed:	Mon Dec 18 19:39:52 2006 UTC (18 years, 6 months ago) by gezelter
Content type:	application/x-tex
File size:	57978 byte(s)
Log Message:	edits