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