ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/multipoleSFPaper/multipoleSFPaper.tex
Revision: 3085
Committed: Mon Dec 18 19:39:52 2006 UTC (18 years, 7 months ago) by gezelter
Content type: application/x-tex
File size: 57978 byte(s)
Log Message:
edits

File Contents

# Content
1 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2 \documentclass[11pt]{article}
3 \usepackage{tabls}
4 %\usepackage{endfloat}
5 \usepackage[tbtags]{amsmath}
6 \usepackage{amsmath,bm}
7 \usepackage{amssymb}
8 \usepackage{mathrsfs}
9 \usepackage{setspace}
10 \usepackage{tabularx}
11 \usepackage{graphicx}
12 \usepackage{booktabs}
13 \usepackage{colortbl}
14 \usepackage[ref]{overcite}
15 \pagestyle{plain}
16 \pagenumbering{arabic}
17 \oddsidemargin 0.0cm \evensidemargin 0.0cm
18 \topmargin -21pt \headsep 10pt
19 \textheight 9.0in \textwidth 6.5in
20 \brokenpenalty=10000
21 \renewcommand{\baselinestretch}{1.2}
22 \renewcommand\citemid{\ } % no comma in optional reference note
23 %\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions
24 %\let\Caption\caption
25 %\renewcommand\caption[1]{%
26 % \Caption[#1]{}%
27 %}
28 %\setlength{\abovecaptionskip}{1.2in}
29 %\setlength{\belowcaptionskip}{1.2in}
30
31 \begin{document}
32
33 \title{Pairwise Alternatives to the Ewald Sum: Applications
34 and Extension to Point Multipoles}
35
36 \author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
37 Department of Chemistry and Biochemistry\\
38 University of Notre Dame\\
39 Notre Dame, Indiana 46556}
40
41 \date{\today}
42
43 \maketitle
44 %\doublespacing
45
46 \begin{abstract}
47 The damped, shifted-force electrostatic potential has been shown to
48 give nearly quantitative agreement with smooth particle mesh Ewald for
49 energy differences between configurations as well as for atomic force
50 and molecular torque vectors. In this paper, we extend this technique
51 to handle interactions between electrostatic multipoles. We also
52 investigate the effects of damped and shifted electrostatics on the
53 static, thermodynamic, and dynamic properties of liquid water and
54 various polymorphs of ice. Additionally, we provide a way of choosing
55 the optimal damping strength for a given cutoff radius that reproduces
56 the static dielectric constant in a variety of water models.
57 \end{abstract}
58
59 \newpage
60
61 %\narrowtext
62
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64 % BODY OF TEXT
65 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66
67 \section{Introduction}
68
69 Over the past several years, there has been increasing interest in
70 pairwise methods for correcting electrostatic interactions in computer
71 simulations of condensed molecular
72 systems.\cite{Wolf99,Zahn02,Kast03,Beck05,Ma05,Fennell06} These
73 techniques were developed from the observations and efforts of Wolf
74 {\it et al.} and their work towards an $\mathcal{O}(N)$ Coulombic
75 sum.\cite{Wolf99} Wolf's method of cutoff neutralization is able to
76 obtain excellent agreement with Madelung energies in ionic
77 crystals.\cite{Wolf99}
78
79 In a recent paper, we showed that simple modifications to Wolf's
80 method could give nearly quantitative agreement with smooth particle
81 mesh Ewald (SPME) for quantities of interest in Monte Carlo
82 (i.e. configurational energy differences) and Molecular Dynamics
83 (i.e. atomic force and molecular torque vectors).\cite{Fennell06} We
84 described the undamped and damped shifted potential (SP) and shifted
85 force (SF) techniques. The potential for the damped form of the SP
86 method, where $\alpha$ is the adjustable damping parameter, is given
87 by
88 \begin{equation}
89 V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}
90 - \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right)
91 \quad r_{ij}\leqslant R_\textrm{c},
92 \label{eq:DSPPot}
93 \end{equation}
94 while the damped form of the SF method is given by
95 \begin{equation}
96 \begin{split}
97 V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
98 \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}
99 - \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\
100 &+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}
101 + \frac{2\alpha}{\pi^{1/2}}
102 \frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
103 \right)\left(r_{ij}-R_\mathrm{c}\right)\ \Biggr{]}
104 \quad r_{ij}\leqslant R_\textrm{c}.
105 \label{eq:DSFPot}
106 \end{split}
107 \end{equation}
108 In these potentials and in all electrostatic quantities that follow,
109 the standard $4 \pi \epsilon_{0}$ has been omitted for clarity.
110
111 The damped SF method is an improvement over the SP method because
112 there is no discontinuity in the forces as particles move out of the
113 cutoff radius ($R_\textrm{c}$). This is accomplished by shifting the
114 forces (and potential) to zero at $R_\textrm{c}$. This is analogous to
115 neutralizing the charge as well as the force effect of the charges
116 within $R_\textrm{c}$.
117
118 To complete the charge neutralization procedure, a self-neutralization
119 term is added to the potential. This term is constant (as long as the
120 charges and cutoff radius do not change), and exists outside the
121 normal pair-loop. It can be thought of as a contribution from a
122 charge opposite in sign, but equal in magnitude, to the central
123 charge, which has been spread out over the surface of the cutoff
124 sphere. This term is calculated via a single loop over all charges in
125 the system. For the undamped case, the self term can be written as
126 \begin{equation}
127 V_\textrm{self} = - \frac{1}{2 R_\textrm{c}} \sum_{i=1}^N q_i^{2},
128 \label{eq:selfTerm}
129 \end{equation}
130 while for the damped case it can be written as
131 \begin{equation}
132 V_\textrm{self} = - \left(\frac{\alpha}{\sqrt{\pi}}
133 + \frac{\textrm{erfc}(\alpha
134 R_\textrm{c})}{2R_\textrm{c}}\right) \sum_{i=1}^N q_i^{2}.
135 \label{eq:dampSelfTerm}
136 \end{equation}
137 The first term within the parentheses of equation
138 (\ref{eq:dampSelfTerm}) is identical to the self term in the Ewald
139 summation, and comes from the utilization of the complimentary error
140 function for electrostatic damping.\cite{deLeeuw80,Wolf99}
141
142 The SF, SP, and Wolf methods operate by neutralizing the total charge
143 contained within the cutoff sphere surrounding each particle. This is
144 accomplished by shifting the potential functions to generate image
145 charges on the surface of the cutoff sphere for each pair interaction
146 computed within $R_\textrm{c}$. The damping function applied to the
147 potential is also an important method for accelerating convergence.
148 In the case of systems involving electrostatic distributions of higher
149 order than point charges, the question remains: How will the shifting
150 and damping need to be modified in order to accommodate point
151 multipoles?
152
153 \section{Electrostatic Damping for Point
154 Multipoles}\label{sec:dampingMultipoles}
155
156 To apply the SF method for systems involving point multipoles, we
157 consider separately the two techniques (shifting and damping) which
158 contribute to the effectiveness of the DSF potential.
159
160 As noted above, shifting the potential and forces is employed to
161 neutralize the total charge contained within each cutoff sphere;
162 however, in a system composed purely of point multipoles, each cutoff
163 sphere is already neutral, so shifting becomes unnecessary. It would
164 still be possible, however, to subtract out the effects of image
165 multipoles. This is essentially what is done for dipoles in the
166 reaction field methods, and the technique could be extended to higher
167 order multipoles.
168
169 In a mixed system of charges and multipoles, the undamped SF potential
170 needs only to shift the force terms between charges and smoothly
171 truncate the multipolar interactions with a switching function. The
172 switching function is required for energy conservation, because a
173 discontinuity will exist in both the potential and forces at
174 $R_\textrm{c}$ in the absence of shifting terms.
175
176 To dampen the SF potential for point multipoles, we need to incorporate
177 the complimentary error function term into the standard forms of the
178 multipolar potentials. We can determine the necessary damping
179 functions by replacing $1/r_{ij}$ with $\mathrm{erfc}(\alpha r_{ij})/r_{ij}$ in the
180 multipole expansion. This procedure quickly becomes quite complex
181 with ``two-center'' systems, like the dipole-dipole potential, and is
182 typically approached using spherical harmonics.\cite{Hirschfelder67} A
183 simpler method for determining damped multipolar interaction
184 potentials arises when we adopt the tensor formalism described by
185 Stone.\cite{Stone02}
186
187 The tensor formalism for electrostatic interactions involves obtaining
188 the multipole interactions from successive gradients of the monopole
189 potential. Thus, tensors of rank zero through two are
190 \begin{equation}
191 T = \frac{1}{r_{ij}},
192 \label{eq:tensorRank1}
193 \end{equation}
194 \begin{equation}
195 T_\alpha = \nabla_\alpha \frac{1}{r_{ij}},
196 \label{eq:tensorRank2}
197 \end{equation}
198 \begin{equation}
199 T_{\alpha\beta} = \nabla_\alpha\nabla_\beta \frac{1}{r_{ij}},
200 \label{eq:tensorRank3}
201 \end{equation}
202 where the form of the first tensor is the charge-charge potential, the
203 second gives the charge-dipole potential, and the third gives the
204 charge-quadrupole and dipole-dipole potentials.\cite{Stone02} Since
205 the force is $-\nabla V$, the forces for each potential come from the
206 next higher tensor.
207
208 As one would do with the multipolar expansion, we can replace $r_{ij}^{-1}$
209 with $\mathrm{erfc}(\alpha r_{ij})/r_{ij}$ to obtain damped forms of the
210 electrostatic potentials. Equation \ref{eq:tensorRank2} generates a
211 damped charge-dipole potential,
212 \begin{equation}
213 V_\textrm{Dcd} = -q_i\frac{\mathbf{r}_{ij}\cdot\boldsymbol{\mu}_j}{r^3_{ij}}
214 c_1(r_{ij}),
215 \label{eq:dChargeDipole}
216 \end{equation}
217 where $c_1(r_{ij})$ is
218 \begin{equation}
219 c_1(r_{ij}) = \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}}
220 + \textrm{erfc}(\alpha r_{ij}).
221 \label{eq:c1Func}
222 \end{equation}
223 Equation \ref{eq:tensorRank3} generates a damped dipole-dipole potential,
224 \begin{equation}
225 V_\textrm{Ddd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
226 (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}}
227 c_2(r_{ij}) -
228 \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}}
229 c_1(r_{ij}),
230 \label{eq:dampDipoleDipole}
231 \end{equation}
232 where
233 \begin{equation}
234 c_2(r_{ij}) = \frac{4\alpha^3r^3_{ij}e^{-\alpha^2r^2_{ij}}}{3\sqrt{\pi}}
235 + \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}}
236 + \textrm{erfc}(\alpha r_{ij}).
237 \end{equation}
238 Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional
239 term. Continuing with higher rank tensors, we can obtain the damping
240 functions for higher multipole potentials and forces. Each subsequent
241 damping function includes one additional term, and we can simplify the
242 procedure for obtaining these terms by writing out the following
243 recurrence relation,
244 \begin{equation}
245 c_n(r_{ij}) = \frac{2^n(\alpha r_{ij})^{2n-1}e^{-\alpha^2r^2_{ij}}}
246 {(2n-1)!!\sqrt{\pi}} + c_{n-1}(r_{ij}),
247 \label{eq:dampingGeneratingFunc}
248 \end{equation}
249 where,
250 \begin{equation}
251 m!! \equiv \left\{ \begin{array}{l@{\quad\quad}l}
252 m\cdot(m-2)\ldots 5\cdot3\cdot1 & m > 0\textrm{ odd} \\
253 m\cdot(m-2)\ldots 6\cdot4\cdot2 & m > 0\textrm{ even} \\
254 1 & m = -1\textrm{ or }0,
255 \end{array}\right.
256 \label{eq:doubleFactorial}
257 \end{equation}
258 and $c_0(r_{ij})$ is erfc$(\alpha r_{ij})$. This generating function
259 is similar in form to those obtained by Smith,\cite{Smith82,Smith98}
260 Toukmaji {\it et al.}\cite{Toukmaji00}, Aguado and
261 Madden,\cite{Aguado03}, and Sagui {\it et al.}\cite{Sagui04} for the
262 application of the Ewald sum to multipoles.
263
264 Returning to the dipole-dipole example, the potential consists of a
265 portion dependent upon $r_{ij}^{-5}$ and another dependent upon
266 $r_{ij}^{-3}$. $c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts
267 respectively. The forces for the damped dipole-dipole interaction, are
268 obtained from the next higher tensor, $T_{\alpha \beta \gamma}$,
269 \begin{equation}
270 \begin{split}
271 F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
272 (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})\mathbf{r}_{ij}}{r^7_{ij}}
273 c_3(r_{ij})\\ &-
274 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})\cdot\boldsymbol{\mu}_j +
275 (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})\cdot\boldsymbol{\mu}_i +
276 \boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}}
277 {r^5_{ij}} c_2(r_{ij}),
278 \end{split}
279 \label{eq:dampDipoleDipoleForces}
280 \end{equation}
281 Using the tensor formalism, we can dampen higher order multipolar
282 interactions using the same effective damping function that we use for
283 charge-charge interactions. This allows us to include multipoles in
284 simulations involving damped electrostatic interactions. In general,
285 if the multipolar potentials are left in $\mathbf{r}_{ij}/r_{ij}$
286 form, instead of reducing them to the more common angular forms which
287 use $\hat{\mathbf{r}}_{ij}$ (or the resultant angles), one may simply replace
288 any $1/r_{ij}^{2n+1}$ dependence with $c_n(r_{ij}) / r_{ij}^{2n+1}$ to
289 obtain the damped version of that multipolar potential.
290
291 As a practical consideration, we note that the evaluation of the
292 complementary error function inside a pair loop can become quite
293 costly. In practice, we pre-compute the $c_n(r)$ functions over a
294 grid of $r$ values and use cubic spline interpolation to obtain
295 estimates of these functions when necessary. Using this procedure,
296 the computational cost of damped electrostatics is only marginally
297 more costly than the undamped case.
298
299 \section{Applications of Damped Shifted-Force Electrostatics}
300
301 Our earlier work on the SF method showed that it can give nearly
302 quantitive agreement with SPME-derived configurational energy
303 differences. The force and torque vectors in identical configurations
304 are also nearly equivalent under the damped SF potential and
305 SPME.\cite{Fennell06} Although these measures bode well for the
306 performance of the SF method in both Monte Carlo and Molecular
307 Dynamics simulations, it would be helpful to have direct comparisons
308 of structural, thermodynamic, and dynamic quantities. To address
309 this, we performed a detailed analysis of a group of simulations
310 involving water models (both point charge and multipolar) under a
311 number of different simulation conditions.
312
313 To provide the most difficult test for the damped SF method, we have
314 chosen a model that has been optimized for use with Ewald sum, and
315 have compared the simulated properties to those computed via Ewald.
316 It is well known that water models parametrized for use with the Ewald
317 sum give calculated properties that are in better agreement with
318 experiment.\cite{vanderSpoel98,Horn04,Rick04} For these reasons, we
319 chose the TIP5P-E water model for our comparisons involving point
320 charges.\cite{Rick04}
321
322 The TIP5P-E water model is a variant of Mahoney and Jorgensen's
323 five-point transferable intermolecular potential (TIP5P) model for
324 water.\cite{Mahoney00} TIP5P was developed to reproduce the density
325 maximum in liquid water near 4$^\circ$C. As with many previous point
326 charge water models (such as ST2, TIP3P, TIP4P, SPC, and SPC/E), TIP5P
327 was parametrized using a simple cutoff with no long-range
328 electrostatic
329 correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87}
330 Without this correction, the pressure term on the central particle
331 from the surroundings is missing. When this correction is included,
332 the system expands to compensate for this added pressure and therefore
333 under-predicts the density of water under standard conditions. In
334 developing TIP5P-E, Rick preserved the geometry and point charge
335 magnitudes in TIP5P and focused on altering the Lennard-Jones
336 parameters to correct the density at 298~K. With the density
337 corrected, he compared common water properties for TIP5P-E using the
338 Ewald sum with TIP5P using a 9~\AA\ cutoff.
339
340 In the following sections, we compare these same properties calculated
341 from TIP5P-E using the Ewald sum with TIP5P-E using the damped SF
342 technique. Our comparisons include the SF technique with a 12~\AA\
343 cutoff and an $\alpha$ = 0.0, 0.1, and 0.2~\AA$^{-1}$, as well as a
344 9~\AA\ cutoff with an $\alpha$ = 0.2~\AA$^{-1}$.
345
346 Moving beyond point-charge electrostatics, the soft sticky dipole
347 (SSD) family of water models is the perfect test case for the
348 point-multipolar extension of damped electrostatics. SSD water models
349 are single point molecules that consist of a ``soft'' Lennard-Jones
350 sphere, a point-dipole, and a tetrahedral function for capturing the
351 hydrogen bonding nature of water - a spherical harmonic term for
352 water-water tetrahedral interactions and a point-quadrupole for
353 interactions with surrounding charges. Detailed descriptions of these
354 models can be found in other
355 studies.\cite{Liu96b,Chandra99,Tan03,Fennell04}
356
357 In deciding which version of the SSD model to use, we need only
358 consider that the SF technique was presented as a pairwise replacement
359 for the Ewald summation. It has been suggested that models
360 parametrized for the Ewald summation (like TIP5P-E) would be
361 appropriate for use with a reaction field and vice versa.\cite{Rick04}
362 Therefore, we decided to test the SSD/RF water model, which was
363 parametrized for use with a reaction field, with the damped
364 electrostatic technique to see how the calculated properties change.
365
366 \subsection{The Density Maximum of TIP5P-E}\label{sec:t5peDensity}
367
368 To compare densities, $NPT$ simulations were performed with the same
369 temperatures as those selected by Rick in his Ewald summation
370 simulations.\cite{Rick04} In order to improve statistics around the
371 density maximum, 3~ns trajectories were accumulated at 0, 12.5, and
372 25$^\circ$C, while 2~ns trajectories were obtained at all other
373 temperatures. The average densities were calculated from the latter
374 three-fourths of each trajectory. Similar to Mahoney and Jorgensen's
375 method for accumulating statistics, these sequences were spliced into
376 200 segments, each providing an average density. These 200 density
377 values were used to calculate the average and standard deviation of
378 the density at each temperature.\cite{Mahoney00}
379
380 \begin{figure}
381 \includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf}
382 \caption{Density versus temperature for the TIP5P-E water model when
383 using the Ewald summation (Ref. \citen{Rick04}) and the SF method with
384 varying cutoff radii and damping coefficients. The pressure term from
385 the image-charge shell is larger than that provided by the
386 reciprocal-space portion of the Ewald summation, leading to slightly
387 lower densities. This effect is more visible with the 9~\AA\ cutoff,
388 where the image charges exert a greater force on the central
389 particle. The representative error bar for the SF methods shows the
390 average one-sigma uncertainty of the density measurement, and this
391 uncertainty is the same for all the SF curves.}
392 \label{fig:t5peDensities}
393 \end{figure}
394 Figure \ref{fig:t5peDensities} shows the densities calculated for
395 TIP5P-E using differing electrostatic corrections overlaid with the
396 experimental values.\cite{CRC80} The densities when using the SF
397 technique are close to, but typically lower than, those calculated
398 using the Ewald summation. These slightly reduced densities indicate
399 that the pressure component from the image charges at R$_\textrm{c}$
400 is larger than that exerted by the reciprocal-space portion of the
401 Ewald summation. Bringing the image charges closer to the central
402 particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than the
403 preferred 12~\AA\ R$_\textrm{c}$) increases the strength of the image
404 charge interactions on the central particle and results in a further
405 reduction of the densities.
406
407 Because the strength of the image charge interactions has a noticeable
408 effect on the density, we would expect the use of electrostatic
409 damping to also play a role in these calculations. Larger values of
410 $\alpha$ weaken the pair-interactions; and since electrostatic damping
411 is distance-dependent, force components from the image charges will be
412 reduced more than those from particles close the the central
413 charge. This effect is visible in figure \ref{fig:t5peDensities} with
414 the damped SF sums showing slightly higher densities than the undamped
415 case; however, it is clear that the choice of cutoff radius plays a
416 much more important role in the resulting densities.
417
418 All of the above density calculations were performed with systems of
419 512 water molecules. Rick observed a system size dependence of the
420 computed densities when using the Ewald summation, most likely due to
421 his tying of the convergence parameter to the box
422 dimensions.\cite{Rick04} For systems of 256 water molecules, the
423 calculated densities were on average 0.002 to 0.003~g/cm$^3$ lower. A
424 system size of 256 molecules would force the use of a shorter
425 R$_\textrm{c}$ when using the SF technique, and this would also lower
426 the densities. Moving to larger systems, as long as the R$_\textrm{c}$
427 remains at a fixed value, we would expect the densities to remain
428 constant.
429
430 \subsection{Liquid Structure of TIP5P-E}\label{sec:t5peLiqStructure}
431
432 The experimentally-determined oxygen-oxygen pair correlation function
433 ($g_\textrm{OO}(r)$) for liquid water has been compared in great
434 detail with predictions of the various common water models, and TIP5P
435 was found to be in better agreement than other rigid, non-polarizable
436 models.\cite{Sorenson00} This excellent agreement with experiment was
437 maintained when Rick developed TIP5P-E.\cite{Rick04} To check whether
438 the choice of using the Ewald summation or the SF technique alters the
439 liquid structure, we calculated this correlation function at 298~K and
440 1~atm for the parameters used in the previous section.
441
442 \begin{figure}
443 \includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf}
444 \caption{The oxygen-oxygen pair correlation functions calculated for
445 TIP5P-E at 298~K and 1~atm while using the Ewald summation
446 (Ref. \citen{Rick04}) and the SF technique with varying
447 parameters. Even with the lower densities obtained using the SF
448 technique, the correlation functions are essentially identical.}
449 \label{fig:t5peGofRs}
450 \end{figure}
451 The pair correlation functions calculated for TIP5P-E while using the
452 SF technique with various parameters are overlaid on the same function
453 obtained while using the Ewald summation in
454 figure~\ref{fig:t5peGofRs}. The differences in density do not appear
455 to have any effect on the liquid structure as the correlation
456 functions are indistinguishable. These results do indicate that
457 $g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic
458 correction.
459
460 \subsection{Thermodynamic Properties of TIP5P-E}\label{sec:t5peThermo}
461
462 In addition to the density and structual features of the liquid, there
463 are a variety of thermodynamic quantities that can be calculated for
464 water and compared directly to experimental values. Some of these
465 additional quantities include the latent heat of vaporization ($\Delta
466 H_\textrm{vap}$), the constant pressure heat capacity ($C_p$), the
467 isothermal compressibility ($\kappa_T$), the thermal expansivity
468 ($\alpha_p$), and the static dielectric constant ($\epsilon$). All of
469 these properties were calculated for TIP5P-E with the Ewald summation,
470 so they provide a good set of reference data for comparisons involving
471 the SF technique.
472
473 The $\Delta H_\textrm{vap}$ is the enthalpy change required to
474 transform one mole of substance from the liquid phase to the gas
475 phase.\cite{Berry00} In molecular simulations, this quantity can be
476 determined via
477 \begin{equation}
478 \begin{split}
479 \Delta H_\textrm{vap} &= H_\textrm{gas} - H_\textrm{liq} \\
480 &= E_\textrm{gas} - E_\textrm{liq}
481 + P(V_\textrm{gas} - V_\textrm{liq}) \\
482 &\approx -\frac{\langle U_\textrm{liq}\rangle}{N}+ RT,
483 \end{split}
484 \label{eq:DeltaHVap}
485 \end{equation}
486 where $E$ is the total energy, $U$ is the potential energy, $P$ is the
487 pressure, $V$ is the volume, $N$ is the number of molecules, $R$ is
488 the gas constant, and $T$ is the temperature.\cite{Jorgensen98b} As
489 seen in the last line of equation (\ref{eq:DeltaHVap}), we can
490 approximate $\Delta H_\textrm{vap}$ by using the ideal gas for the gas
491 state. This allows us to cancel the kinetic energy terms, leaving only
492 the $U_\textrm{liq}$ term. Additionally, the $PV$ term for the gas is
493 several orders of magnitude larger than that of the liquid, so we can
494 neglect the liquid $PV$ term.
495
496 The remaining thermodynamic properties can all be calculated from
497 fluctuations of the enthalpy, volume, and system dipole
498 moment.\cite{Allen87} $C_p$ can be calculated from fluctuations in the
499 enthalpy in constant pressure simulations via
500 \begin{equation}
501 \begin{split}
502 C_p = \left(\frac{\partial H}{\partial T}\right)_{N,P}
503 = \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2},
504 \end{split}
505 \label{eq:Cp}
506 \end{equation}
507 where $k_B$ is Boltzmann's constant. $\kappa_T$ can be calculated via
508 \begin{equation}
509 \begin{split}
510 \kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T}
511 = \frac{({\langle V^2\rangle}_{NPT} - {\langle V\rangle}^{2}_{NPT})}
512 {k_BT\langle V\rangle_{NPT}},
513 \end{split}
514 \label{eq:kappa}
515 \end{equation}
516 and $\alpha_p$ can be calculated via
517 \begin{equation}
518 \begin{split}
519 \alpha_p = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P}
520 = \frac{(\langle VH\rangle_{NPT}
521 - \langle V\rangle_{NPT}\langle H\rangle_{NPT})}
522 {k_BT^2\langle V\rangle_{NPT}}.
523 \end{split}
524 \label{eq:alpha}
525 \end{equation}
526 Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can
527 be calculated for systems of non-polarizable substances via
528 \begin{equation}
529 \epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV},
530 \label{eq:staticDielectric}
531 \end{equation}
532 where $\epsilon_0$ is the permittivity of free space and $\langle
533 M^2\rangle$ is the fluctuation of the system dipole
534 moment.\cite{Allen87} The numerator in the fractional term in equation
535 (\ref{eq:staticDielectric}) is the fluctuation of the simulation-box
536 dipole moment, identical to the quantity calculated in the
537 finite-system Kirkwood $g$ factor ($G_k$):
538 \begin{equation}
539 G_k = \frac{\langle M^2\rangle}{N\mu^2},
540 \label{eq:KirkwoodFactor}
541 \end{equation}
542 where $\mu$ is the dipole moment of a single molecule of the
543 homogeneous system.\cite{Neumann83,Neumann84,Neumann85} The
544 fluctuation term in both equation (\ref{eq:staticDielectric}) and
545 (\ref{eq:KirkwoodFactor}) is calculated as follows,
546 \begin{equation}
547 \begin{split}
548 \langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle
549 - \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\
550 &= \langle M_x^2+M_y^2+M_z^2\rangle
551 - (\langle M_x\rangle^2 + \langle M_x\rangle^2
552 + \langle M_x\rangle^2).
553 \end{split}
554 \label{eq:fluctBoxDipole}
555 \end{equation}
556 This fluctuation term can be accumulated during the simulation;
557 however, it converges rather slowly, thus requiring multi-nanosecond
558 simulation times.\cite{Horn04} In the case of tin-foil boundary
559 conditions, the dielectric/surface term of the Ewald summation is
560 equal to zero. Since the SF method also lacks this
561 dielectric/surface term, equation (\ref{eq:staticDielectric}) is still
562 valid for determining static dielectric constants.
563
564 All of the above properties were calculated from the same trajectories
565 used to determine the densities in section \ref{sec:t5peDensity}
566 except for the static dielectric constants. The $\epsilon$ values were
567 accumulated from 2~ns $NVE$ ensemble trajectories with system densities
568 fixed at the average values from the $NPT$ simulations at each of the
569 temperatures. The resulting values are displayed in figure
570 \ref{fig:t5peThermo}.
571 \begin{figure}
572 \centering
573 \includegraphics[width=5.8in]{./figures/t5peThermo.pdf}
574 \caption{Thermodynamic properties for TIP5P-E using the Ewald summation
575 and the SF techniques along with the experimental values. Units
576 for the properties are kcal mol$^{-1}$ for $\Delta H_\textrm{vap}$,
577 cal mol$^{-1}$ K$^{-1}$ for $C_p$, 10$^6$ atm$^{-1}$ for $\kappa_T$,
578 and 10$^5$ K$^{-1}$ for $\alpha_p$. Ewald summation results are from
579 reference \citen{Rick04}. Experimental values for $\Delta
580 H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference
581 \citen{Kell75}. Experimental values for $C_p$ are from reference
582 \citen{Wagner02}. Experimental values for $\epsilon$ are from reference
583 \citen{Malmberg56}.}
584 \label{fig:t5peThermo}
585 \end{figure}
586
587 For all of the properties computed, the trends with temperature
588 obtained when using the Ewald summation are reproduced with the SF
589 technique. One noticeable difference between the properties calculated
590 using the two methods are the lower $\Delta H_\textrm{vap}$ values
591 when using SF. This is to be expected due to the direct weakening of
592 the electrostatic interaction through forced neutralization. This
593 results in an increase of the intermolecular potential producing lower
594 values from equation (\ref{eq:DeltaHVap}). The slopes of these values
595 with temperature are similar to that seen using the Ewald summation;
596 however, they are both steeper than the experimental trend, indirectly
597 resulting in the inflated $C_p$ values at all temperatures.
598
599 Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ values
600 all overlap within error. As indicated for the $\Delta H_\textrm{vap}$
601 and $C_p$ results, the deviations between experiment and simulation in
602 this region are not the fault of the electrostatic summation methods
603 but are due to the geometry and parameters of the TIP5P class of water
604 models. Like most rigid, non-polarizable, point-charge water models,
605 the density decreases with temperature at a much faster rate than
606 experiment (see figure \ref{fig:t5peDensities}). This reduced density
607 leads to the inflated compressibility and expansivity values at higher
608 temperatures seen here in figure \ref{fig:t5peThermo}. Incorporation
609 of polarizability and many-body effects are required in order for
610 water models to overcome differences between simulation-based and
611 experimentally determined densities at these higher
612 temperatures.\cite{Laasonen93,Donchev06}
613
614 At temperatures below the freezing point for experimental water, the
615 differences between SF and the Ewald summation results are more
616 apparent. The larger $C_p$ and lower $\alpha_p$ values in this region
617 indicate a more pronounced transition in the supercooled regime,
618 particularly in the case of SF without damping. This points to the
619 onset of a more frustrated or glassy behavior for the undamped and
620 weakly-damped SF simulations of TIP5P-E at temperatures below 250~K
621 than is seen from the Ewald sum at these temperatures. Undamped SF
622 electrostatics has a stronger contribution from nearby charges.
623 Damping these local interactions or using a reciprocal-space method
624 makes the water less sensitive to ordering on a shorter length scale.
625 We can recover nearly quantitative agreement with the Ewald results by
626 increasing the damping constant.
627
628 The final thermodynamic property displayed in figure
629 \ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
630 between the Ewald and SF methods (and with experiment). It is known
631 that the dielectric constant is dependent upon and is quite sensitive
632 to the imposed boundary conditions.\cite{Neumann80,Neumann83} This is
633 readily apparent in the converged $\epsilon$ values accumulated for
634 the SF simulations. Lack of a damping function results in dielectric
635 constants significantly smaller than those obtained using the Ewald
636 sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the
637 agreement considerably.
638
639 It should be noted that the choice of the ``Ewald coefficient''
640 ($\kappa$) and real-space cutoff values also have a significant effect
641 on the calculated static dielectric constant when using the Ewald
642 summation. In general, aqueous systems with larger screening
643 parameters will report larger values for the dielectric constant, the
644 same behavior we see here with SF. Although the received wisdom
645 (CITATION) on the Ewald parameter is that large enough Fourier-space
646 grids can make up for the effects of overdamping the real-space
647 portion of the Ewald sum, this wisdom does not appear to be borne out
648 by actual numerical results. To illustrate this point, we have
649 computed the static dielectric constants from 216-molecule SPC water
650 simulations (9 \AA\ cutoff radius, 3 ns data collection times, PME
651 with fourth-order spline to smooth the reciprocal space portion of the
652 sum). We carried out these simulations with a range of different
653 Ewald coefficients. For each value of the Ewald coefficient, we also
654 varied the number of Fourier vectors used to compute the $k$-space
655 sum. At the largest Fourier-grid (48x48x48 Fourier-vectors), this
656 corresponds to a 0.39 \AA\ grid spacing. The results for these
657 calculations are shown in figure \ref{fig:ewaldDielectric}.
658
659 \begin{figure}
660 \includegraphics[width=\linewidth]{./figures/ewaldDielectric.pdf}
661 \caption{The dielectric constant of the SPC water model using the Ewald
662 summation as a function of ({\it left}) the Ewald coefficient
663 ($\kappa$) and ({\it right}) the $k$-Space grid detail. All
664 calculations were performed with a 216 water molecule system using a
665 fixed cutoff radius of 9 \AA\ at 300 K.}
666 \label{fig:ewaldDielectric}
667 \end{figure}
668
669 It is clear that the choice of the Ewald coefficient actually does
670 have a dramatic effect on the static dielectric of water, and that
671 increasing the Fourier grid to 3 times the typical value does not
672 bring these dielectric constants back into agreement with each other.
673 It should be noted that this behavior would be troubling in {\it any}
674 method for calculating electrostatic interactions. Ideally, a method
675 should give fixed values for important properties like the static
676 dielectric and should be insensitive to the choice of convergence
677 parameters like the Ewald coefficient and $\alpha$. From what we have
678 observed, neither the Ewald summation nor the real-space shifted
679 methods appear to be ``black boxes'' at this point. Both require
680 careful consideration of the choice of parameters on the property
681 being studied.
682
683 \subsection{Optimal Damping Coefficients for Damped
684 Electrostatics}\label{sec:t5peDielectric}
685
686 In the previous section, we observed that the choice of damping
687 coefficient plays a major role in the calculated dielectric constant
688 for the SF method. Similar damping parameter behavior was observed in
689 the long-time correlated motions of the NaCl crystal.\cite{Fennell06}
690 The static dielectric constant is calculated from the long-time
691 fluctuations of the system's accumulated dipole moment
692 (Eq. (\ref{eq:staticDielectric})), so it is quite sensitive to the
693 choice of damping parameter. Since $\alpha$ is an adjustable
694 parameter, it would be best to choose optimal damping constants such
695 that any arbitrary choice of cutoff radius will properly capture the
696 dielectric behavior of the liquid.
697
698 In order to find these optimal values, we mapped out the static
699 dielectric constant as a function of both the damping parameter and
700 cutoff radius for TIP5P-E and for a point-dipolar water model
701 (SSD/RF). To calculate the static dielectric constant, we performed
702 5~ns $NPT$ calculations on systems of 512 water molecules and averaged
703 over the converged region (typically the latter 2.5~ns) of equation
704 (\ref{eq:staticDielectric}). The selected cutoff radii ranged from 9,
705 10, 11, and 12~\AA , and the damping parameter values ranged from 0.1
706 to 0.35~\AA$^{-1}$.
707
708 \begin{table}
709 \centering
710 \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
711 cutoff radius $R_\textrm{c}$. The color scale ranges from blue (10) to red (100).}
712 \vspace{6pt}
713 \begin{tabular}{ lccccccccc }
714 \toprule
715 \toprule
716 & \multicolumn{4}{c}{TIP5P-E} & & \multicolumn{4}{c}{SSD/RF} \\
717 & \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} & & \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\
718 \cmidrule(lr){2-5} \cmidrule(lr){7-10}
719 $\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 & & 9 & 10 & 11 & 12 \\
720 \midrule
721 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 \\
722 & \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} \\
723 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 \\
724 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 \\
725 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 \\
726 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 \\
727 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 \\
728 & \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} \\
729 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 \\
730 & \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} \\
731 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 \\
732 \bottomrule
733 \end{tabular}
734 \label{tab:DielectricMap}
735 \end{table}
736
737 The results of these calculations are displayed in table
738 \ref{tab:DielectricMap}. The dielectric constants for both models
739 decrease with increasing cutoff radii ($R_\textrm{c}$) and increase
740 with increasing damping ($\alpha$). Another point to note is that
741 choosing $\alpha$ and $R_\textrm{c}$ identical to those used with the
742 Ewald summation results in the same calculated dielectric constant. As
743 an example, in the paper outlining the development of TIP5P-E, the
744 real-space cutoff and Ewald coefficient were tethered to the system
745 size, and for a 512 molecule system are approximately 12~\AA\ and
746 0.25~\AA$^{-1}$ respectively.\cite{Rick04} These parameters resulted
747 in a dielectric constant of 92$\pm$14, while with SF these parameters
748 give a dielectric constant of 90.8$\pm$0.9. Another example comes from
749 the TIP4P-Ew paper where $\alpha$ and $R_\textrm{c}$ were chosen to be
750 9.5~\AA\ and 0.35~\AA$^{-1}$, and these parameters resulted in a
751 dielectric constant equal to 63$\pm$1.\cite{Horn04} Calculations using
752 SF with these parameters and this water model give a dielectric
753 constant of 61$\pm$1. Since the dielectric constant is dependent on
754 $\alpha$ and $R_\textrm{c}$ with the SF technique, it might be
755 interesting to investigate the dependence of the static dielectric
756 constant on the choice of convergence parameters ($R_\textrm{c}$ and
757 $\kappa$) utilized in most implementations of the Ewald sum.
758
759 It is also apparent from this table that electrostatic damping has a
760 more pronounced effect on the dipolar interactions of SSD/RF than the
761 monopolar interactions of TIP5P-E. The dielectric constant covers a
762 much wider range and has a steeper slope with increasing damping
763 parameter.
764
765 Although it is tempting to choose damping parameters equivalent to the
766 Ewald examples to obtain quantitative agreement, the results of our
767 previous work indicate that values this high are destructive to both
768 the energetics and dynamics.\cite{Fennell06} Ideally, $\alpha$ should
769 not exceed 0.3~\AA$^{-1}$ for any of the cutoff values in this
770 range.
771
772 It is possible to use a simple energetic criterion for choosing an
773 optimal value of $\alpha$ given a cutoff radius ($R_\textrm{c}$) and
774 an energetic tolerance. This is identical to how the Ewald
775 coefficient ($\kappa$) is chosen in many of the common simulation
776 packages (GROMACS,\cite{Gromacs} AMBER,\cite{AMBER} and
777 TINKER~\cite{TINKER}). In the case of the damped method described in
778 this work, the value of the tolerance,
779 \begin{equation}
780 \textrm{tolerance} = \frac{V_\textrm{damped}(R_\textrm{c})}{V_\textrm{coulomb}(R_\textrm{c})}
781 = \textrm{erfc}(\alpha R_\textrm{c}),
782 \end{equation}
783 which does best at reproducing the static dielectric constant in a
784 variety of water models is $3.1~\times~10^{-4}$. Using this
785 tolerance, good choices of $\alpha$ for cutoff radii of 9 and 12
786 \AA\ are 0.2834 and 0.2125 \AA$^{-1}$, respectively.
787
788 As a final note on optimal damping parameters, aside from a slight
789 lowering of $\Delta H_\textrm{vap}$, using optimal $\alpha$ values
790 does not alter any of the other thermodynamic properties.
791
792 \subsection{Dynamic Properties of TIP5P-E}\label{sec:t5peDynamics}
793
794 To look at the dynamic properties of TIP5P-E when using the SF method,
795 200~ps $NVE$ simulations were performed for each temperature at the
796 average density obtained from the $NPT$ simulations. $R_\textrm{c}$
797 values of 9 and 12~\AA\ and the ideal $\alpha$ values determined above
798 (0.2834 and 0.2125~\AA$^{-1}$) were used for the damped
799 electrostatics. The self-diffusion constants (D) were calculated from
800 linear fits to the long-time portion of the mean square displacement
801 ($\langle r^{2}(t) \rangle$).\cite{Allen87}
802
803 In addition to translational diffusion, orientational relaxation times
804 were calculated for comparisons with the Ewald simulations and with
805 experiments. These values were determined from the same 200~ps $NVE$
806 trajectories used for translational diffusion by calculating the
807 orientational time correlation function,
808 \begin{equation}
809 C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\gamma(t)
810 \cdot\hat{\mathbf{u}}_i^\gamma(0)\right]\right\rangle,
811 \label{eq:OrientCorr}
812 \end{equation}
813 where $P_l$ is the Legendre polynomial of order $l$ and
814 $\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along
815 axis $\gamma$. The body-fixed reference frame used for our
816 orientational correlation functions has the $z$-axis running along the
817 HOH bisector, and the $y$-axis connecting the two hydrogen atoms.
818 $C_l^y$ is therefore calculated from the time evolution of a vector of
819 unit length pointing between the two hydrogen atoms.
820
821 From the orientation autocorrelation functions, we can obtain time
822 constants for rotational relaxation. The relatively short time
823 portions (between 1 and 3~ps for water) of these curves can be fit to
824 an exponential decay to obtain these constants, and they are directly
825 comparable to water orientational relaxation times from nuclear
826 magnetic resonance (NMR). The relaxation constant obtained from
827 $C_2^y(t)$ is of particular interest because it describes the
828 relaxation of the principle axis connecting the hydrogen atoms. Thus,
829 $C_2^y(t)$ can be compared to the intermolecular portion of the
830 dipole-dipole relaxation from a proton NMR signal and should provide
831 the best estimate of the NMR relaxation time constant.\cite{Impey82}
832
833 \begin{figure}
834 \centering
835 \includegraphics[width=5.8in]{./figures/t5peDynamics.pdf}
836 \caption{Diffusion constants ({\it upper}) and reorientational time
837 constants ({\it lower}) for TIP5P-E using the Ewald sum and SF
838 technique compared with experiment. Data at temperatures less than
839 0$^\circ$C were not included in the $\tau_2^y$ plot to allow for
840 easier comparisons in the more relevant temperature regime.}
841 \label{fig:t5peDynamics}
842 \end{figure}
843 Results for the diffusion constants and orientational relaxation times
844 are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
845 apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
846 the Ewald sum are reproduced with the SF technique. The enhanced
847 diffusion (relative to experiment) at high temperatures are again the
848 product of the lower simulated densities and do not provide any
849 special insight into differences between the electrostatic summation
850 techniques. Though not apparent in this figure, SF values at the
851 lowest temperature are approximately twice as slow as $D$ values
852 obtained using the Ewald sum. These values support the observation
853 from section \ref{sec:t5peThermo} that the SF simulations result in a
854 slightly more viscous supercooled region than is obtained using the
855 Ewald sum.
856
857 The $\tau_2^y$ results in the lower frame of figure
858 \ref{fig:t5peDynamics} show a much greater difference between the SF
859 results and the Ewald results. At all temperatures shown, TIP5P-E
860 relaxes faster than experiment with the Ewald sum while tracking
861 experiment fairly well when using the SF technique (independent of the
862 choice of damping constant). There are several possible reasons for
863 this deviation between techniques. The Ewald results were calculated
864 using shorter trajectories (10~ps) than the SF results (200~ps).
865 Calculation of these SF values from 10~ps trajectories (with
866 subsequently lower accuracy) showed a 0.4~ps drop in $\tau_2^y$,
867 placing the result more in line with that obtained using the Ewald
868 sum. Recomputing correlation times to meet a lower statistical
869 standard is counter-productive, however. Assuming the Ewald results
870 are not entirely the product of poor statistics, differences in
871 techniques to integrate the orientational motion could also play a
872 role. {\sc shake} is the most commonly used technique for
873 approximating rigid-body orientational motion,\cite{Ryckaert77}
874 whereas in {\sc oopse}, we maintain and integrate the entire rotation
875 matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake}
876 is an iterative constraint technique, if the convergence tolerances
877 are raised for increased performance, error will accumulate in the
878 orientational motion. Finally, the Ewald results were calculated using
879 the $NVT$ ensemble, while the $NVE$ ensemble was used for SF
880 calculations. The motion due to the extended variable (the thermostat)
881 will always alter the dynamics, resulting in differences between $NVT$
882 and $NVE$ results. These differences are increasingly noticeable as
883 the time constant for the thermostat decreases.
884
885 \subsection{Comparison of Reaction Field and Damped Electrostatics for
886 SSD/RF}
887
888 SSD/RF was parametrized for use with a reaction field, which is a
889 common and relatively inexpensive way of handling long-range
890 electrostatic corrections in dipolar systems.\cite{Onsager36}
891 Although there is no reason to expect that damped electrostatics will
892 behave in a similar fashion to the reaction field, it is well known
893 that model that are parametrized for use with Ewald do better than
894 unoptimized models under the influence of a reaction
895 field.\cite{Rick04} We compared a number of properties of SSD/RF that
896 had previously been computed using a reaction field with those same
897 values under damped electrostatics.
898
899 The properties shown in table \ref{tab:dampedSSDRF} show that using
900 damped electrostatics can result in even better agreement with
901 experiment than is obtained via reaction field. The average density
902 shows a modest increase when using damped electrostatics in place of
903 the reaction field. This comes about because we neglect the pressure
904 effect due to the surroundings outside of the cutoff, instead relying
905 on screening effects to neutralize electrostatic interactions at long
906 distances. The $C_p$ also shows a slight increase, indicating greater
907 fluctuation in the enthalpy at constant pressure. The only other
908 differences between the damped electrostatics and the reaction field
909 results are the dipole reorientational time constants, $\tau_1$ and
910 $\tau_2$. When using damped electrostatics, the water molecules relax
911 more quickly and exhibit behavior closer to the experimental
912 values. These results indicate that since there is no need to specify
913 an external dielectric constant with the damped electrostatics, it is
914 almost certainly a better choice for dipolar simulations than the
915 reaction field method. Additionally, by using damped electrostatics
916 instead of reaction field, SSD/RF can be used effectively with mixed
917 charge / dipolar systems, such as dissolved ions, dissolved organic
918 molecules, or even proteins.
919
920 \begin{table}
921 \caption{Properties of SSD/RF when using reaction field or damped
922 electrostatics. Simulations were carried out at 298~K, 1~atm, and
923 with 512 molecules.}
924 \footnotesize
925 \centering
926 \begin{tabular}{ llccc }
927 \toprule
928 \toprule
929 & & Reaction Field (Ref. \citen{Fennell04}) & Damped Electrostatics &
930 Experiment [Ref.] \\
931 & & $\epsilon = 80$ & $R_\textrm{c} = 12$\AA ; $\alpha = 0.2125$~\AA$^{-1}$ & \\
932 \midrule
933 $\rho$ & (g cm$^{-3}$) & 0.997(1) & 1.004(4) & 0.997 [\citen{CRC80}]\\
934 $C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 23.8(2) & 27(1) & 18.005 [\citen{Wagner02}] \\
935 $D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.32(6) & 2.33(2) & 2.299 [\citen{Mills73}]\\
936 $n_C$ & & 4.4 & 4.2 & 4.7 [\citen{Hura00}]\\
937 $n_H$ & & 3.7 & 3.7 & 3.5 [\citen{Soper86}]\\
938 $\tau_1$ & (ps) & 7.2(4) & 5.86(8) & 5.7 [\citen{Eisenberg69}]\\
939 $\tau_2$ & (ps) & 3.2(2) & 2.45(7) & 2.3 [\citen{Krynicki66}]\\
940 \bottomrule
941 \end{tabular}
942 \label{tab:dampedSSDRF}
943 \end{table}
944
945 \subsection{Predictions of Ice Polymorph Stability}
946
947 In an earlier paper, we performed a series of free energy calculations
948 on several ice polymorphs which are stable or metastable at low
949 pressures, one of which (Ice-$i$) we observed in spontaneous
950 crystallizations of an early version of the SSD/RF water
951 model.\cite{Fennell05} In this study, a distinct dependence of the
952 computed free energies on the cutoff radius and electrostatic
953 summation method was observed. Since the SF technique can be used as
954 a simple and efficient replacement for the Ewald summation, our final
955 test of this method is to see if it is capable of addressing the
956 spurious stability of the Ice-$i$ phases with the various common water
957 models. To this end, we have performed thermodynamic integrations of
958 all of the previously discussed ice polymorphs using the SF technique
959 with a cutoff radius of 12~\AA\ and an $\alpha$ of 0.2125~\AA . These
960 calculations were performed on TIP5P-E and TIP4P-Ew (variants of the
961 TIP5P and TIP4P models optimized for the Ewald summation) as well as
962 SPC/E and SSD/RF.
963
964 \begin{table}
965 \centering
966 \caption{Helmholtz free energies of ice polymorphs at 1~atm and 200~K
967 using the damped SF electrostatic correction method with a
968 variety of water models.}
969 \begin{tabular}{ lccccc }
970 \toprule
971 \toprule
972 Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\
973 \cmidrule(lr){2-6}
974 & \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\
975 \midrule
976 TIP5P-E & -11.98(4) & -11.96(4) & -11.87(3) & - & -11.95(3) \\
977 TIP4P-Ew & -13.11(3) & -13.09(3) & -12.97(3) & - & -12.98(3) \\
978 SPC/E & -12.99(3) & -13.00(3) & -13.03(3) & - & -12.99(3) \\
979 SSD/RF & -11.83(3) & -11.66(4) & -12.32(3) & -12.39(3) & - \\
980 \bottomrule
981 \end{tabular}
982 \label{tab:dampedFreeEnergy}
983 \end{table}
984 The results of these calculations in table \ref{tab:dampedFreeEnergy}
985 show similar behavior to the Ewald results in the previous
986 study.\cite{Fennell05} The Helmholtz free energies of the ice
987 polymorphs with SSD/RF order in the same fashion, with Ice-$i$ having
988 the lowest free energy; however, the Ice-$i$ and ice B free energies
989 are quite a bit closer (nearly isoenergetic). The SPC/E results show
990 the different polymorphs to be nearly isoenergetic. This is the same
991 behavior observed using an Ewald correction.\cite{Fennell05} Ice B has
992 the lowest Helmholtz free energy for SPC/E; however, all the polymorph
993 results overlap within the error estimates.
994
995 The most interesting results from these calculations come from the
996 more expensive TIP4P-Ew and TIP5P-E results. Both of these models were
997 optimized for use with an electrostatic correction and are
998 geometrically arranged to mimic water using drastically different
999 charge distributions. In TIP5P-E, the primary location for the
1000 negative charge in the molecule is assigned to the lone-pairs of the
1001 oxygen, while TIP4P-Ew places the negative charge near the
1002 center-of-mass along the H-O-H bisector. There is some debate as to
1003 which is the proper choice for the negative charge location, and this
1004 has in part led to a six-site water model that balances both of these
1005 options.\cite{Vega05,Nada03} The limited results in table
1006 \ref{tab:dampedFreeEnergy} support the results of Vega {\it et al.},
1007 which indicate the TIP4P charge location geometry performs better
1008 under some circumstances.\cite{Vega05} With the TIP4P-Ew water model,
1009 the experimentally observed polymorph (ice I$_\textrm{h}$) is the
1010 preferred form with ice I$_\textrm{c}$ slightly higher in energy,
1011 though overlapping within error. Additionally, the spurious ice B and
1012 Ice-$i^\prime$ structures are destabilized relative to these
1013 polymorphs. TIP5P-E shows similar behavior to SPC/E, where there is no
1014 real free energy distinction between the various polymorphs. While ice
1015 B is close in free energy to the other polymorphs, these results fail
1016 to support the findings of other researchers indicating the preferred
1017 form of TIP5P at 1~atm is a structure similar to ice
1018 B.\cite{Yamada02,Vega05,Abascal05} It should be noted that we were
1019 looking at TIP5P-E rather than TIP5P, and the differences in the
1020 Lennard-Jones parameters could cause this discrepancy. Overall, these
1021 results indicate that TIP4P-Ew is a better mimic of the solid forms of
1022 water than some of the other models.
1023
1024 \section{Conclusions}
1025
1026 This investigation of pairwise electrostatic summation techniques
1027 shows that there is a viable and computationally efficient alternative
1028 to the Ewald summation. The SF method (equation (\ref{eq:DSFPot}))
1029 has proven itself capable of reproducing structural, thermodynamic,
1030 and dynamic quantities that are nearly quantitative matches to results
1031 from far more expensive methods. Additionally, we have now extended
1032 the damping formalism to electrostatic multipoles, so the damped SF
1033 potential can be used in systems that contain mixtures of charges and
1034 point multipoles.
1035
1036 We have also provided a simple prescription for choosing optimal
1037 damping parameters given a choice of cutoff radius. The damping
1038 parameters were chosen to obtain a static dielectric constant as close
1039 as possible to the experimental value, which should be useful for
1040 simulating the electrostatic screening properties of liquid water
1041 accurately. The formula for optimal damping was the same for a
1042 complicated multipoint model as it was for a simple point-dipolar
1043 model, and is identical to energetic tolerance methods commonly used
1044 to choose the Ewald coefficient.
1045
1046 As in all purely pairwise cutoff methods, the damped SF method is
1047 expected to scale approximately {\it linearly} with system size, and
1048 is easily parallelizable. This should result in substantial
1049 reductions in the computational cost of performing large simulations.
1050 With the proper use of pre-computation and spline interpolation, the
1051 damped SF method is essentially the same cost as a simple real-space
1052 cutoff.
1053
1054 We are not suggesting that there is any flaw with the Ewald sum; in
1055 fact, it is the standard by which the damped SF method has been
1056 judged. However, these results provide further evidence that in the
1057 typical simulations performed today, the Ewald summation may no longer
1058 be required to obtain the level of accuracy most researchers have come
1059 to expect.
1060
1061 \section{Acknowledgments}
1062 Support for this project was provided by the National Science
1063 Foundation under grant CHE-0134881. Computation time was provided by
1064 the Notre Dame Center for Research Computing. The authors would like
1065 to thank Steve Corcelli and Ed Maginn for helpful discussions and
1066 comments.
1067
1068 \newpage
1069
1070 \bibliographystyle{achemso}
1071 \bibliography{multipoleSFPaper}
1072
1073
1074 \end{document}