1 |
chrisfen |
3064 |
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
2 |
gezelter |
3066 |
\documentclass[11pt]{article} |
3 |
chrisfen |
3064 |
\usepackage{tabls} |
4 |
gezelter |
3068 |
%\usepackage{endfloat} |
5 |
chrisfen |
3064 |
\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 |
gezelter |
3068 |
%\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 |
chrisfen |
3064 |
|
31 |
|
|
\begin{document} |
32 |
|
|
|
33 |
gezelter |
3066 |
\title{Pairwise Alternatives to the Ewald Sum: Applications |
34 |
|
|
and Extension to Point Multipoles} |
35 |
chrisfen |
3064 |
|
36 |
gezelter |
3068 |
\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 |
chrisfen |
3064 |
Notre Dame, Indiana 46556} |
40 |
|
|
|
41 |
|
|
\date{\today} |
42 |
|
|
|
43 |
|
|
\maketitle |
44 |
|
|
%\doublespacing |
45 |
|
|
|
46 |
|
|
\begin{abstract} |
47 |
gezelter |
3066 |
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. We provide a way of choosing the optimal |
55 |
|
|
damping strength for a given cutoff radius that reproduces the static |
56 |
|
|
dielectric constant in a variety of water models. |
57 |
chrisfen |
3064 |
\end{abstract} |
58 |
|
|
|
59 |
gezelter |
3068 |
\newpage |
60 |
|
|
|
61 |
chrisfen |
3064 |
%\narrowtext |
62 |
|
|
|
63 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
64 |
|
|
% BODY OF TEXT |
65 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
66 |
|
|
|
67 |
|
|
\section{Introduction} |
68 |
|
|
|
69 |
gezelter |
3066 |
Over the past several years, there has been increasing interest in |
70 |
|
|
pairwise methods for correcting electrostatic interactions in computer |
71 |
chrisfen |
3064 |
simulations of condensed molecular |
72 |
gezelter |
3066 |
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 |
chrisfen |
3064 |
|
79 |
gezelter |
3066 |
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 |
gezelter |
3067 |
force (SF) techniques. The potential for the damped form of the SP |
86 |
|
|
method is given by |
87 |
chrisfen |
3064 |
\begin{equation} |
88 |
gezelter |
3067 |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r} |
89 |
|
|
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) |
90 |
|
|
\quad r\leqslant R_\textrm{c}, |
91 |
|
|
\label{eq:DSPPot} |
92 |
|
|
\end{equation} |
93 |
|
|
while the damped form of the SF method is given by |
94 |
|
|
\begin{equation} |
95 |
chrisfen |
3064 |
\begin{split} |
96 |
|
|
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
97 |
|
|
\frac{\mathrm{erfc}\left(\alpha r\right)}{r} |
98 |
|
|
- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ |
99 |
|
|
&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} |
100 |
|
|
+ \frac{2\alpha}{\pi^{1/2}} |
101 |
|
|
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} |
102 |
|
|
\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]} |
103 |
|
|
\quad r\leqslant R_\textrm{c}. |
104 |
|
|
\label{eq:DSFPot} |
105 |
|
|
\end{split} |
106 |
|
|
\end{equation} |
107 |
gezelter |
3067 |
(In these potentials and in all electrostatic quantities that follow, |
108 |
gezelter |
3066 |
the standard $4 \pi \epsilon_{0}$ has been omitted for clarity.) |
109 |
chrisfen |
3064 |
|
110 |
gezelter |
3066 |
The damped SF method is an improvement over the SP method because |
111 |
|
|
there is no discontinuity in the forces as particles move out of the |
112 |
|
|
cutoff radius ($R_\textrm{c}$). This is accomplished by shifting the |
113 |
|
|
forces (and potential) to zero at $R_\textrm{c}$. This is analogous to |
114 |
|
|
neutralizing the charge as well as the force effect of the charges |
115 |
|
|
within $R_\textrm{c}$. |
116 |
|
|
|
117 |
|
|
To complete the charge neutralization procedure, a self-neutralization |
118 |
gezelter |
3067 |
term is included in the potential. This term is constant (as long as |
119 |
|
|
the charges and cutoff radius do not change), and exists outside the |
120 |
|
|
normal pair-loop. It can be thought of as a contribution from a |
121 |
|
|
charge opposite in sign, but equal in magnitude, to the central |
122 |
gezelter |
3066 |
charge, but which has been spread out over the surface of the cutoff |
123 |
|
|
sphere. This term is calculated via a single loop over all charges in |
124 |
|
|
the system. For the undamped case, the self term can be written as |
125 |
chrisfen |
3064 |
\begin{equation} |
126 |
gezelter |
3066 |
V_\textrm{self} = \frac{1}{2 R_\textrm{c}} \sum_{i=1}^N q_i^{2}, |
127 |
chrisfen |
3064 |
\label{eq:selfTerm} |
128 |
|
|
\end{equation} |
129 |
|
|
while for the damped case it can be written as |
130 |
|
|
\begin{equation} |
131 |
gezelter |
3066 |
V_\textrm{self} = \left(\frac{\alpha}{\sqrt{\pi}} |
132 |
|
|
+ \frac{\textrm{erfc}(\alpha |
133 |
|
|
R_\textrm{c})}{2R_\textrm{c}}\right) \sum_{i=1}^N q_i^{2}. |
134 |
chrisfen |
3064 |
\label{eq:dampSelfTerm} |
135 |
|
|
\end{equation} |
136 |
|
|
The first term within the parentheses of equation |
137 |
gezelter |
3066 |
(\ref{eq:dampSelfTerm}) is identical to the self term in the Ewald |
138 |
chrisfen |
3064 |
summation, and comes from the utilization of the complimentary error |
139 |
|
|
function for electrostatic damping.\cite{deLeeuw80,Wolf99} |
140 |
|
|
|
141 |
gezelter |
3066 |
The SF, SP, and Wolf methods operate by neutralizing the total charge |
142 |
|
|
contained within the cutoff sphere surrounding each particle. This is |
143 |
|
|
accomplished by creating image charges on the surface of the cutoff |
144 |
|
|
sphere for each pair interaction computed within the sphere. The |
145 |
|
|
damping function applied to the potential is also an important method |
146 |
|
|
for accelerating convergence. In the case of systems involving |
147 |
|
|
electrostatic distributions of higher order than point charges, the |
148 |
|
|
question remains: How will the shifting and damping need to be |
149 |
|
|
modified in order to accommodate point multipoles? |
150 |
chrisfen |
3064 |
|
151 |
gezelter |
3066 |
\section{Electrostatic Damping for Point |
152 |
|
|
Multipoles}\label{sec:dampingMultipoles} |
153 |
chrisfen |
3064 |
|
154 |
gezelter |
3066 |
To apply the SF method for systems involving point multipoles, we |
155 |
|
|
consider separately the two techniques (shifting and damping) which |
156 |
|
|
contribute to the effectiveness of the DSF potential. |
157 |
chrisfen |
3064 |
|
158 |
gezelter |
3066 |
As noted above, shifting the potential and forces is employed to |
159 |
|
|
neutralize the total charge contained within each cutoff sphere; |
160 |
|
|
however, in a system composed purely of point multipoles, each cutoff |
161 |
|
|
sphere is already neutral, so shifting becomes unnecessary. |
162 |
|
|
|
163 |
|
|
In a mixed system of charges and multipoles, the undamped SF potential |
164 |
|
|
needs only to shift the force terms between charges and smoothly |
165 |
|
|
truncate the multipolar interactions with a switching function. The |
166 |
|
|
switching function is required for energy conservation, because a |
167 |
|
|
discontinuity will exist in both the potential and forces at |
168 |
|
|
$R_\textrm{c}$ in the absence of shifting terms. |
169 |
|
|
|
170 |
|
|
To damp the SF potential for point multipoles, we need to incorporate |
171 |
|
|
the complimentary error function term into the standard forms of the |
172 |
|
|
multipolar potentials. We can determine the necessary damping |
173 |
|
|
functions by replacing $1/r$ with $\mathrm{erfc}(\alpha r)/r$ in the |
174 |
|
|
multipole expansion. This procedure quickly becomes quite complex |
175 |
|
|
with ``two-center'' systems, like the dipole-dipole potential, and is |
176 |
|
|
typically approached using spherical harmonics.\cite{Hirschfelder67} A |
177 |
|
|
simpler method for determining damped multipolar interaction |
178 |
|
|
potentials arises when we adopt the tensor formalism described by |
179 |
|
|
Stone.\cite{Stone02} |
180 |
|
|
|
181 |
|
|
The tensor formalism for electrostatic interactions involves obtaining |
182 |
|
|
the multipole interactions from successive gradients of the monopole |
183 |
|
|
potential. Thus, tensors of rank zero through two are |
184 |
chrisfen |
3064 |
\begin{equation} |
185 |
gezelter |
3066 |
T = \frac{1}{r_{ij}}, |
186 |
|
|
\label{eq:tensorRank1} |
187 |
chrisfen |
3064 |
\end{equation} |
188 |
|
|
\begin{equation} |
189 |
gezelter |
3066 |
T_\alpha = \nabla_\alpha \frac{1}{r_{ij}}, |
190 |
|
|
\label{eq:tensorRank2} |
191 |
chrisfen |
3064 |
\end{equation} |
192 |
gezelter |
3066 |
\begin{equation} |
193 |
|
|
T_{\alpha\beta} = \nabla_\alpha\nabla_\beta \frac{1}{r_{ij}}, |
194 |
|
|
\label{eq:tensorRank3} |
195 |
|
|
\end{equation} |
196 |
|
|
where the form of the first tensor is the charge-charge potential, the |
197 |
|
|
second gives the charge-dipole potential, and the third gives the |
198 |
|
|
charge-quadrupole and dipole-dipole potentials.\cite{Stone02} Since |
199 |
|
|
the force is $-\nabla V$, the forces for each potential come from the |
200 |
|
|
next higher tensor. |
201 |
chrisfen |
3064 |
|
202 |
gezelter |
3066 |
As one would do with the multipolar expansion, we can replace $r^{-1}$ |
203 |
gezelter |
3067 |
with $\mathrm{erfc}(\alpha r)/r$ to obtain damped forms of the |
204 |
gezelter |
3066 |
electrostatic potentials. Equation \ref{eq:tensorRank2} generates a |
205 |
|
|
damped charge-dipole potential, |
206 |
chrisfen |
3064 |
\begin{equation} |
207 |
gezelter |
3068 |
V_\textrm{Dcd} = -q_i\frac{\mathbf{r}_{ij}\cdot\boldsymbol{\mu}_j}{r^3_{ij}} |
208 |
|
|
c_1(r_{ij}), |
209 |
chrisfen |
3064 |
\label{eq:dChargeDipole} |
210 |
|
|
\end{equation} |
211 |
|
|
where $c_1(r_{ij})$ is |
212 |
|
|
\begin{equation} |
213 |
|
|
c_1(r_{ij}) = \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}} |
214 |
|
|
+ \textrm{erfc}(\alpha r_{ij}). |
215 |
|
|
\label{eq:c1Func} |
216 |
|
|
\end{equation} |
217 |
|
|
|
218 |
gezelter |
3066 |
Equation \ref{eq:tensorRank3} generates a damped dipole-dipole potential, |
219 |
chrisfen |
3064 |
\begin{equation} |
220 |
gezelter |
3066 |
V_\textrm{Ddd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}) |
221 |
|
|
(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}} |
222 |
|
|
c_2(r_{ij}) - |
223 |
|
|
\frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}} |
224 |
|
|
c_1(r_{ij}), |
225 |
|
|
\label{eq:dampDipoleDipole} |
226 |
chrisfen |
3064 |
\end{equation} |
227 |
gezelter |
3066 |
where |
228 |
chrisfen |
3064 |
\begin{equation} |
229 |
|
|
c_2(r_{ij}) = \frac{4\alpha^3r^3_{ij}e^{-\alpha^2r^2_{ij}}}{3\sqrt{\pi}} |
230 |
|
|
+ \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}} |
231 |
|
|
+ \textrm{erfc}(\alpha r_{ij}). |
232 |
|
|
\end{equation} |
233 |
|
|
Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional |
234 |
|
|
term. Continuing with higher rank tensors, we can obtain the damping |
235 |
|
|
functions for higher multipole potentials and forces. Each subsequent |
236 |
|
|
damping function includes one additional term, and we can simplify the |
237 |
|
|
procedure for obtaining these terms by writing out the following |
238 |
|
|
generating function, |
239 |
|
|
\begin{equation} |
240 |
|
|
c_n(r_{ij}) = \frac{2^n(\alpha r_{ij})^{2n-1}e^{-\alpha^2r^2_{ij}}} |
241 |
|
|
{(2n-1)!!\sqrt{\pi}} + c_{n-1}(r_{ij}), |
242 |
|
|
\label{eq:dampingGeneratingFunc} |
243 |
|
|
\end{equation} |
244 |
|
|
where, |
245 |
|
|
\begin{equation} |
246 |
|
|
m!! \equiv \left\{ \begin{array}{l@{\quad\quad}l} |
247 |
|
|
m\cdot(m-2)\ldots 5\cdot3\cdot1 & m > 0\textrm{ odd} \\ |
248 |
|
|
m\cdot(m-2)\ldots 6\cdot4\cdot2 & m > 0\textrm{ even} \\ |
249 |
|
|
1 & m = -1\textrm{ or }0, |
250 |
|
|
\end{array}\right. |
251 |
|
|
\label{eq:doubleFactorial} |
252 |
|
|
\end{equation} |
253 |
|
|
and $c_0(r_{ij})$ is erfc$(\alpha r_{ij})$. This generating function |
254 |
gezelter |
3066 |
is similar in form to those obtained by Smith and Aguado and Madden |
255 |
|
|
for the application of the Ewald sum to |
256 |
chrisfen |
3064 |
multipoles.\cite{Smith82,Smith98,Aguado03} |
257 |
|
|
|
258 |
|
|
Returning to the dipole-dipole example, the potential consists of a |
259 |
|
|
portion dependent upon $r^{-5}$ and another dependent upon |
260 |
gezelter |
3066 |
$r^{-3}$. |
261 |
chrisfen |
3064 |
$c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts |
262 |
gezelter |
3066 |
respectively. The forces for the damped dipole-dipole interaction, are |
263 |
|
|
obtained from the next higher tensor, $T_{\alpha \beta \gamma}$, |
264 |
chrisfen |
3064 |
\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\hat{\boldsymbol{\mu}}_j + |
270 |
|
|
\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}\cdot\hat{\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 |
gezelter |
3066 |
Using the tensor formalism, we can dampen higher order multipolar |
277 |
|
|
interactions using the same effective damping function that we use for |
278 |
gezelter |
3068 |
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{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 |
chrisfen |
3064 |
|
286 |
gezelter |
3068 |
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 |
chrisfen |
3064 |
|
294 |
gezelter |
3068 |
\section{Applications of Damped Shifted-Force Electrostatics} |
295 |
gezelter |
3067 |
|
296 |
gezelter |
3066 |
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 |
chrisfen |
3064 |
|
308 |
gezelter |
3066 |
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 |
chrisfen |
3064 |
|
317 |
gezelter |
3066 |
The soft sticky dipole (SSD) family of water models is the perfect |
318 |
|
|
test case for the point-multipolar extension of damped electrostatics. |
319 |
|
|
SSD water models are single point molecules that consist of a ``soft'' |
320 |
|
|
Lennard-Jones sphere, a point-dipole, and a tetrahedral function for |
321 |
|
|
capturing the hydrogen bonding nature of water - a spherical harmonic |
322 |
|
|
term for water-water tetrahedral interactions and a point-quadrupole |
323 |
|
|
for interactions with surrounding charges. Detailed descriptions of |
324 |
|
|
these models can be found in other |
325 |
|
|
studies.\cite{Liu96b,Chandra99,Tan03,Fennell04} |
326 |
|
|
|
327 |
|
|
In deciding which version of the SSD model to use, we need only |
328 |
|
|
consider that the SF technique was presented as a pairwise replacement |
329 |
|
|
for the Ewald summation. It has been suggested that models |
330 |
|
|
parametrized for the Ewald summation (like TIP5P-E) would be |
331 |
|
|
appropriate for use with a reaction field and vice versa.\cite{Rick04} |
332 |
|
|
Therefore, we decided to test the SSD/RF water model, which was |
333 |
|
|
parametrized for use with a reaction field, with the damped |
334 |
|
|
electrostatic technique to see how the calculated properties change. |
335 |
|
|
|
336 |
chrisfen |
3064 |
The TIP5P-E water model is a variant of Mahoney and Jorgensen's |
337 |
|
|
five-point transferable intermolecular potential (TIP5P) model for |
338 |
|
|
water.\cite{Mahoney00} TIP5P was developed to reproduce the density |
339 |
gezelter |
3066 |
maximum in liquid water near 4$^\circ$C. As with many previous point |
340 |
|
|
charge water models (such as ST2, TIP3P, TIP4P, SPC, and SPC/E), TIP5P |
341 |
|
|
was parametrized using a simple cutoff with no long-range |
342 |
|
|
electrostatic |
343 |
chrisfen |
3064 |
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87} |
344 |
|
|
Without this correction, the pressure term on the central particle |
345 |
|
|
from the surroundings is missing. When this correction is included, |
346 |
|
|
systems of these particles expand to compensate for this added |
347 |
|
|
pressure term and under-predict the density of water under standard |
348 |
gezelter |
3066 |
conditions. In developing TIP5P-E, Rick preserved the geometry and |
349 |
|
|
point charge magnitudes in TIP5P and focused on altering the |
350 |
|
|
Lennard-Jones parameters to correct the density at 298~K. With the |
351 |
chrisfen |
3064 |
density corrected, he compared common water properties for TIP5P-E |
352 |
|
|
using the Ewald sum with TIP5P using a 9~\AA\ cutoff. |
353 |
|
|
|
354 |
gezelter |
3066 |
In the following sections, we compare these same properties calculated |
355 |
|
|
from TIP5P-E using the Ewald sum with TIP5P-E using the damped SF |
356 |
|
|
technique. Our comparisons include the SF technique with a 12~\AA\ |
357 |
|
|
cutoff and an $\alpha$ = 0.0, 0.1, and 0.2~\AA$^{-1}$, as well as a |
358 |
|
|
9~\AA\ cutoff with an $\alpha$ = 0.2~\AA$^{-1}$. |
359 |
chrisfen |
3064 |
|
360 |
gezelter |
3066 |
\subsection{The Density Maximum of TIP5P-E}\label{sec:t5peDensity} |
361 |
chrisfen |
3064 |
|
362 |
|
|
To compare densities, $NPT$ simulations were performed with the same |
363 |
|
|
temperatures as those selected by Rick in his Ewald summation |
364 |
|
|
simulations.\cite{Rick04} In order to improve statistics around the |
365 |
|
|
density maximum, 3~ns trajectories were accumulated at 0, 12.5, and |
366 |
|
|
25$^\circ$C, while 2~ns trajectories were obtained at all other |
367 |
|
|
temperatures. The average densities were calculated from the later |
368 |
|
|
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
369 |
|
|
method for accumulating statistics, these sequences were spliced into |
370 |
|
|
200 segments, each providing an average density. These 200 density |
371 |
|
|
values were used to calculate the average and standard deviation of |
372 |
|
|
the density at each temperature.\cite{Mahoney00} |
373 |
|
|
|
374 |
|
|
\begin{figure} |
375 |
|
|
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
376 |
|
|
\caption{Density versus temperature for the TIP5P-E water model when |
377 |
gezelter |
3066 |
using the Ewald summation (Ref. \citen{Rick04}) and the SF method with |
378 |
|
|
varying cutoff radii and damping coefficients. The pressure term from |
379 |
|
|
the image-charge shell is larger than that provided by the |
380 |
|
|
reciprocal-space portion of the Ewald summation, leading to slightly |
381 |
|
|
lower densities. This effect is more visible with the 9~\AA\ cutoff, |
382 |
|
|
where the image charges exert a greater force on the central |
383 |
|
|
particle. The error bars for the SF methods show the average one-sigma |
384 |
|
|
uncertainty of the density measurement, and this uncertainty is the |
385 |
|
|
same for all the SF curves.} |
386 |
chrisfen |
3064 |
\label{fig:t5peDensities} |
387 |
|
|
\end{figure} |
388 |
|
|
Figure \ref{fig:t5peDensities} shows the densities calculated for |
389 |
gezelter |
3066 |
TIP5P-E using differing electrostatic corrections overlaid with the |
390 |
|
|
experimental values.\cite{CRC80} The densities when using the SF |
391 |
|
|
technique are close to, but typically lower than, those calculated |
392 |
chrisfen |
3064 |
using the Ewald summation. These slightly reduced densities indicate |
393 |
|
|
that the pressure component from the image charges at R$_\textrm{c}$ |
394 |
|
|
is larger than that exerted by the reciprocal-space portion of the |
395 |
gezelter |
3066 |
Ewald summation. Bringing the image charges closer to the central |
396 |
chrisfen |
3064 |
particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than the |
397 |
|
|
preferred 12~\AA\ R$_\textrm{c}$) increases the strength of the image |
398 |
|
|
charge interactions on the central particle and results in a further |
399 |
|
|
reduction of the densities. |
400 |
|
|
|
401 |
|
|
Because the strength of the image charge interactions has a noticeable |
402 |
|
|
effect on the density, we would expect the use of electrostatic |
403 |
|
|
damping to also play a role in these calculations. Larger values of |
404 |
|
|
$\alpha$ weaken the pair-interactions; and since electrostatic damping |
405 |
|
|
is distance-dependent, force components from the image charges will be |
406 |
|
|
reduced more than those from particles close the the central |
407 |
|
|
charge. This effect is visible in figure \ref{fig:t5peDensities} with |
408 |
gezelter |
3066 |
the damped SF sums showing slightly higher densities; however, it is |
409 |
|
|
clear that the choice of cutoff radius plays a much more important |
410 |
|
|
role in the resulting densities. |
411 |
chrisfen |
3064 |
|
412 |
gezelter |
3066 |
All of the above density calculations were performed with systems of |
413 |
|
|
512 water molecules. Rick observed a system size dependence of the |
414 |
|
|
computed densities when using the Ewald summation, most likely due to |
415 |
|
|
his tying of the convergence parameter to the box |
416 |
chrisfen |
3064 |
dimensions.\cite{Rick04} For systems of 256 water molecules, the |
417 |
|
|
calculated densities were on average 0.002 to 0.003~g/cm$^3$ lower. A |
418 |
|
|
system size of 256 molecules would force the use of a shorter |
419 |
gezelter |
3066 |
R$_\textrm{c}$ when using the SF technique, and this would also lower |
420 |
|
|
the densities. Moving to larger systems, as long as the R$_\textrm{c}$ |
421 |
|
|
remains at a fixed value, we would expect the densities to remain |
422 |
|
|
constant. |
423 |
chrisfen |
3064 |
|
424 |
gezelter |
3066 |
\subsection{Liquid Structure of TIP5P-E}\label{sec:t5peLiqStructure} |
425 |
chrisfen |
3064 |
|
426 |
gezelter |
3068 |
The experimentally-determined oxygen-oxygen pair correlation function |
427 |
|
|
($g_\textrm{OO}(r)$) for liquid water has been compared in great |
428 |
|
|
detail with predictions of the various common water models, and TIP5P |
429 |
|
|
was found to be in better agreement than other rigid, non-polarizable |
430 |
|
|
models.\cite{Sorenson00} This excellent agreement with experiment was |
431 |
|
|
maintained when Rick developed TIP5P-E.\cite{Rick04} To check whether |
432 |
|
|
the choice of using the Ewald summation or the SF technique alters the |
433 |
|
|
liquid structure, we calculated this correlation function at 298~K and |
434 |
|
|
1~atm for the parameters used in the previous section. |
435 |
chrisfen |
3064 |
|
436 |
|
|
\begin{figure} |
437 |
|
|
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
438 |
gezelter |
3068 |
\caption{The oxygen-oxygen pair correlation functions calculated for |
439 |
|
|
TIP5P-E at 298~K and 1atm while using the Ewald summation |
440 |
|
|
(Ref. \cite{Rick04}) and the SF technique with varying |
441 |
|
|
parameters. Even with the lower densities obtained using the SF |
442 |
|
|
technique, the correlation functions are essentially identical.} |
443 |
chrisfen |
3064 |
\label{fig:t5peGofRs} |
444 |
|
|
\end{figure} |
445 |
gezelter |
3068 |
The pair correlation functions calculated for TIP5P-E while using the |
446 |
|
|
SF technique with various parameters are overlaid on the same function |
447 |
|
|
obtained while using the Ewald summation in |
448 |
chrisfen |
3064 |
figure~\ref{fig:t5peGofRs}. The differences in density do not appear |
449 |
gezelter |
3068 |
to have any effect on the liquid structure as the correlation |
450 |
|
|
functions are indistinguishable. These results do indicate that |
451 |
chrisfen |
3064 |
$g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic |
452 |
|
|
correction. |
453 |
|
|
|
454 |
gezelter |
3066 |
\subsection{Thermodynamic Properties of TIP5P-E}\label{sec:t5peThermo} |
455 |
chrisfen |
3064 |
|
456 |
gezelter |
3066 |
In addition to the density and structual features of the liquid, there |
457 |
|
|
are a variety of thermodynamic quantities that can be calculated for |
458 |
|
|
water and compared directly to experimental values. Some of these |
459 |
|
|
additional quantities include the latent heat of vaporization ($\Delta |
460 |
|
|
H_\textrm{vap}$), the constant pressure heat capacity ($C_p$), the |
461 |
|
|
isothermal compressibility ($\kappa_T$), the thermal expansivity |
462 |
|
|
($\alpha_p$), and the static dielectric constant ($\epsilon$). All of |
463 |
|
|
these properties were calculated for TIP5P-E with the Ewald summation, |
464 |
gezelter |
3068 |
so they provide a good set of reference data for comparisons involving |
465 |
|
|
the SF technique. |
466 |
chrisfen |
3064 |
|
467 |
|
|
The $\Delta H_\textrm{vap}$ is the enthalpy change required to |
468 |
|
|
transform one mole of substance from the liquid phase to the gas |
469 |
|
|
phase.\cite{Berry00} In molecular simulations, this quantity can be |
470 |
|
|
determined via |
471 |
|
|
\begin{equation} |
472 |
|
|
\begin{split} |
473 |
gezelter |
3066 |
\Delta H_\textrm{vap} &= H_\textrm{gas} - H_\textrm{liq} \\ |
474 |
|
|
&= E_\textrm{gas} - E_\textrm{liq} |
475 |
|
|
+ P(V_\textrm{gas} - V_\textrm{liq}) \\ |
476 |
|
|
&\approx -\frac{\langle U_\textrm{liq}\rangle}{N}+ RT, |
477 |
chrisfen |
3064 |
\end{split} |
478 |
|
|
\label{eq:DeltaHVap} |
479 |
|
|
\end{equation} |
480 |
gezelter |
3066 |
where $E$ is the total energy, $U$ is the potential energy, $P$ is the |
481 |
chrisfen |
3064 |
pressure, $V$ is the volume, $N$ is the number of molecules, $R$ is |
482 |
|
|
the gas constant, and $T$ is the temperature.\cite{Jorgensen98b} As |
483 |
|
|
seen in the last line of equation (\ref{eq:DeltaHVap}), we can |
484 |
|
|
approximate $\Delta H_\textrm{vap}$ by using the ideal gas for the gas |
485 |
|
|
state. This allows us to cancel the kinetic energy terms, leaving only |
486 |
gezelter |
3068 |
the $U_\textrm{liq}$ term. Additionally, the $PV$ term for the gas is |
487 |
chrisfen |
3064 |
several orders of magnitude larger than that of the liquid, so we can |
488 |
gezelter |
3068 |
neglect the liquid $PV$ term. |
489 |
chrisfen |
3064 |
|
490 |
|
|
The remaining thermodynamic properties can all be calculated from |
491 |
|
|
fluctuations of the enthalpy, volume, and system dipole |
492 |
|
|
moment.\cite{Allen87} $C_p$ can be calculated from fluctuations in the |
493 |
|
|
enthalpy in constant pressure simulations via |
494 |
|
|
\begin{equation} |
495 |
|
|
\begin{split} |
496 |
gezelter |
3066 |
C_p = \left(\frac{\partial H}{\partial T}\right)_{N,P} |
497 |
chrisfen |
3064 |
= \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2}, |
498 |
|
|
\end{split} |
499 |
|
|
\label{eq:Cp} |
500 |
|
|
\end{equation} |
501 |
|
|
where $k_B$ is Boltzmann's constant. $\kappa_T$ can be calculated via |
502 |
|
|
\begin{equation} |
503 |
|
|
\begin{split} |
504 |
|
|
\kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T} |
505 |
gezelter |
3068 |
= \frac{(\langle V^2\rangle_{NPT} - \langle V\rangle^{2}_{NPT})} |
506 |
gezelter |
3066 |
{k_BT\langle V\rangle_{NPT}}, |
507 |
chrisfen |
3064 |
\end{split} |
508 |
|
|
\label{eq:kappa} |
509 |
|
|
\end{equation} |
510 |
|
|
and $\alpha_p$ can be calculated via |
511 |
|
|
\begin{equation} |
512 |
|
|
\begin{split} |
513 |
|
|
\alpha_p = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P} |
514 |
gezelter |
3066 |
= \frac{(\langle VH\rangle_{NPT} |
515 |
|
|
- \langle V\rangle_{NPT}\langle H\rangle_{NPT})} |
516 |
|
|
{k_BT^2\langle V\rangle_{NPT}}. |
517 |
chrisfen |
3064 |
\end{split} |
518 |
|
|
\label{eq:alpha} |
519 |
|
|
\end{equation} |
520 |
|
|
Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can |
521 |
|
|
be calculated for systems of non-polarizable substances via |
522 |
|
|
\begin{equation} |
523 |
|
|
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV}, |
524 |
|
|
\label{eq:staticDielectric} |
525 |
|
|
\end{equation} |
526 |
|
|
where $\epsilon_0$ is the permittivity of free space and $\langle |
527 |
|
|
M^2\rangle$ is the fluctuation of the system dipole |
528 |
|
|
moment.\cite{Allen87} The numerator in the fractional term in equation |
529 |
|
|
(\ref{eq:staticDielectric}) is the fluctuation of the simulation-box |
530 |
|
|
dipole moment, identical to the quantity calculated in the |
531 |
|
|
finite-system Kirkwood $g$ factor ($G_k$): |
532 |
|
|
\begin{equation} |
533 |
|
|
G_k = \frac{\langle M^2\rangle}{N\mu^2}, |
534 |
|
|
\label{eq:KirkwoodFactor} |
535 |
|
|
\end{equation} |
536 |
|
|
where $\mu$ is the dipole moment of a single molecule of the |
537 |
|
|
homogeneous system.\cite{Neumann83,Neumann84,Neumann85} The |
538 |
|
|
fluctuation term in both equation (\ref{eq:staticDielectric}) and |
539 |
|
|
\ref{eq:KirkwoodFactor} is calculated as follows, |
540 |
|
|
\begin{equation} |
541 |
|
|
\begin{split} |
542 |
|
|
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle |
543 |
|
|
- \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\ |
544 |
|
|
&= \langle M_x^2+M_y^2+M_z^2\rangle |
545 |
|
|
- (\langle M_x\rangle^2 + \langle M_x\rangle^2 |
546 |
|
|
+ \langle M_x\rangle^2). |
547 |
|
|
\end{split} |
548 |
|
|
\label{eq:fluctBoxDipole} |
549 |
|
|
\end{equation} |
550 |
|
|
This fluctuation term can be accumulated during the simulation; |
551 |
|
|
however, it converges rather slowly, thus requiring multi-nanosecond |
552 |
|
|
simulation times.\cite{Horn04} In the case of tin-foil boundary |
553 |
|
|
conditions, the dielectric/surface term of the Ewald summation is |
554 |
gezelter |
3066 |
equal to zero. Since the SF method also lacks this |
555 |
chrisfen |
3064 |
dielectric/surface term, equation (\ref{eq:staticDielectric}) is still |
556 |
|
|
valid for determining static dielectric constants. |
557 |
|
|
|
558 |
|
|
All of the above properties were calculated from the same trajectories |
559 |
|
|
used to determine the densities in section \ref{sec:t5peDensity} |
560 |
|
|
except for the static dielectric constants. The $\epsilon$ values were |
561 |
|
|
accumulated from 2~ns $NVE$ ensemble trajectories with system densities |
562 |
|
|
fixed at the average values from the $NPT$ simulations at each of the |
563 |
|
|
temperatures. The resulting values are displayed in figure |
564 |
|
|
\ref{fig:t5peThermo}. |
565 |
|
|
\begin{figure} |
566 |
|
|
\centering |
567 |
gezelter |
3067 |
\includegraphics[width=5.8in]{./figures/t5peThermo.pdf} |
568 |
chrisfen |
3064 |
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation |
569 |
gezelter |
3066 |
and the SF techniques along with the experimental values. Units |
570 |
chrisfen |
3064 |
for the properties are kcal mol$^{-1}$ for $\Delta H_\textrm{vap}$, |
571 |
|
|
cal mol$^{-1}$ K$^{-1}$ for $C_p$, 10$^6$ atm$^{-1}$ for $\kappa_T$, |
572 |
|
|
and 10$^5$ K$^{-1}$ for $\alpha_p$. Ewald summation results are from |
573 |
|
|
reference \cite{Rick04}. Experimental values for $\Delta |
574 |
|
|
H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference |
575 |
|
|
\cite{Kell75}. Experimental values for $C_p$ are from reference |
576 |
|
|
\cite{Wagner02}. Experimental values for $\epsilon$ are from reference |
577 |
|
|
\cite{Malmberg56}.} |
578 |
|
|
\label{fig:t5peThermo} |
579 |
|
|
\end{figure} |
580 |
|
|
|
581 |
gezelter |
3066 |
For all of the properties computed, the trends with temperature |
582 |
|
|
obtained when using the Ewald summation are reproduced with the SF |
583 |
|
|
technique. One noticeable difference between the properties calculated |
584 |
|
|
using the two methods are the lower $\Delta H_\textrm{vap}$ values |
585 |
|
|
when using SF. This is to be expected due to the direct weakening of |
586 |
|
|
the electrostatic interaction through forced neutralization. This |
587 |
|
|
results in an increase of the intermolecular potential producing lower |
588 |
|
|
values from equation (\ref{eq:DeltaHVap}). The slopes of these values |
589 |
|
|
with temperature are similar to that seen using the Ewald summation; |
590 |
|
|
however, they are both steeper than the experimental trend, indirectly |
591 |
|
|
resulting in the inflated $C_p$ values at all temperatures. |
592 |
chrisfen |
3064 |
|
593 |
|
|
Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$ values |
594 |
|
|
all overlap within error. As indicated for the $\Delta H_\textrm{vap}$ |
595 |
gezelter |
3066 |
and $C_p$ results, the deviations between experiment and simulation in |
596 |
|
|
this region are not the fault of the electrostatic summation methods |
597 |
|
|
but are due to the geometry and parameters of the TIP5P class of water |
598 |
|
|
models. Like most rigid, non-polarizable, point-charge water models, |
599 |
|
|
the density decreases with temperature at a much faster rate than |
600 |
|
|
experiment (see figure \ref{fig:t5peDensities}). This reduced density |
601 |
|
|
leads to the inflated compressibility and expansivity values at higher |
602 |
|
|
temperatures seen here in figure \ref{fig:t5peThermo}. Incorporation |
603 |
|
|
of polarizability and many-body effects are required in order for |
604 |
|
|
water models to overcome differences between simulation-based and |
605 |
|
|
experimentally determined densities at these higher |
606 |
chrisfen |
3064 |
temperatures.\cite{Laasonen93,Donchev06} |
607 |
|
|
|
608 |
|
|
At temperatures below the freezing point for experimental water, the |
609 |
gezelter |
3066 |
differences between SF and the Ewald summation results are more |
610 |
chrisfen |
3064 |
apparent. The larger $C_p$ and lower $\alpha_p$ values in this region |
611 |
|
|
indicate a more pronounced transition in the supercooled regime, |
612 |
gezelter |
3068 |
particularly in the case of SF without damping. This points to the |
613 |
|
|
onset of a more frustrated or glassy behavior for the undamped and |
614 |
|
|
weakly-damped SF simulations of TIP5P-E at temperatures below 250~K |
615 |
|
|
than is seen from the Ewald sum. Undamped SF electrostatics has a |
616 |
|
|
stronger contribution from nearby charges. Damping these local |
617 |
|
|
interactions or using a reciprocal-space method makes the water less |
618 |
|
|
sensitive to ordering on a shorter length scale. We can recover |
619 |
|
|
nearly quantitative agreement with the Ewald results by increasing the |
620 |
|
|
damping constant. |
621 |
chrisfen |
3064 |
|
622 |
|
|
The final thermodynamic property displayed in figure |
623 |
|
|
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy |
624 |
gezelter |
3066 |
between the Ewald and SF methods (and with experiment). It is known |
625 |
|
|
that the dielectric constant is dependent upon and is quite sensitive |
626 |
|
|
to the imposed boundary conditions.\cite{Neumann80,Neumann83} This is |
627 |
|
|
readily apparent in the converged $\epsilon$ values accumulated for |
628 |
|
|
the SF simulations. Lack of a damping function results in dielectric |
629 |
chrisfen |
3064 |
constants significantly smaller than those obtained using the Ewald |
630 |
|
|
sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the |
631 |
|
|
agreement considerably. It should be noted that the choice of the |
632 |
gezelter |
3066 |
``Ewald coefficient'' ($\kappa$) and real-space cutoff values also |
633 |
|
|
have a significant effect on the calculated static dielectric constant |
634 |
|
|
when using the Ewald summation. In the simulations of TIP5P-E with the |
635 |
|
|
Ewald sum, this screening parameter was tethered to the simulation box |
636 |
|
|
size (as was the $R_\textrm{c}$).\cite{Rick04} In general, systems |
637 |
|
|
with larger screening parameters reported larger dielectric constant |
638 |
|
|
values, the same behavior we see here with {\sc sf}; however, the |
639 |
|
|
choice of cutoff radius also plays an important role. |
640 |
chrisfen |
3064 |
|
641 |
gezelter |
3068 |
\subsection{Optimal Damping Coefficients for Damped |
642 |
|
|
Electrostatics}\label{sec:t5peDielectric} |
643 |
chrisfen |
3064 |
|
644 |
gezelter |
3066 |
In the previous section, we observed that the choice of damping |
645 |
|
|
coefficient plays a major role in the calculated dielectric constant |
646 |
|
|
for the SF method. Similar damping parameter behavior was observed in |
647 |
|
|
the long-time correlated motions of the NaCl crystal.\cite{Fennell06} |
648 |
|
|
The static dielectric constant is calculated from the long-time |
649 |
|
|
fluctuations of the system's accumulated dipole moment |
650 |
|
|
(Eq. (\ref{eq:staticDielectric})), so it is quite sensitive to the |
651 |
|
|
choice of damping parameter. Since $\alpha$ is an adjustable |
652 |
|
|
parameter, it would be best to choose optimal damping constants such |
653 |
|
|
that any arbitrary choice of cutoff radius will properly capture the |
654 |
|
|
dielectric behavior of the liquid. |
655 |
chrisfen |
3064 |
|
656 |
|
|
In order to find these optimal values, we mapped out the static |
657 |
|
|
dielectric constant as a function of both the damping parameter and |
658 |
gezelter |
3066 |
cutoff radius for TIP5P-E and for a point-dipolar water model |
659 |
|
|
(SSD/RF). To calculate the static dielectric constant, we performed |
660 |
|
|
5~ns $NPT$ calculations on systems of 512 water molecules and averaged |
661 |
|
|
over the converged region (typically the later 2.5~ns) of equation |
662 |
|
|
(\ref{eq:staticDielectric}). The selected cutoff radii ranged from 9, |
663 |
|
|
10, 11, and 12~\AA , and the damping parameter values ranged from 0.1 |
664 |
|
|
to 0.35~\AA$^{-1}$. |
665 |
chrisfen |
3064 |
|
666 |
|
|
\begin{table} |
667 |
|
|
\centering |
668 |
gezelter |
3066 |
\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 |
669 |
|
|
cutoff radius $R_\textrm{c}$. The color scale ranges from blue (10) to red (100).} |
670 |
chrisfen |
3064 |
\vspace{6pt} |
671 |
gezelter |
3066 |
\begin{tabular}{ lccccccccc } |
672 |
chrisfen |
3064 |
\toprule |
673 |
|
|
\toprule |
674 |
gezelter |
3066 |
& \multicolumn{4}{c}{TIP5P-E} & & \multicolumn{4}{c}{SSD/RF} \\ |
675 |
|
|
& \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} & & \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\ |
676 |
|
|
\cmidrule(lr){2-5} \cmidrule(lr){7-10} |
677 |
|
|
$\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 & & 9 & 10 & 11 & 12 \\ |
678 |
chrisfen |
3064 |
\midrule |
679 |
gezelter |
3066 |
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 \\ |
680 |
|
|
& \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} \\ |
681 |
|
|
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 \\ |
682 |
|
|
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 \\ |
683 |
|
|
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 \\ |
684 |
|
|
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 \\ |
685 |
|
|
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 \\ |
686 |
|
|
& \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} \\ |
687 |
|
|
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 \\ |
688 |
|
|
& \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} \\ |
689 |
|
|
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 \\ |
690 |
chrisfen |
3064 |
\bottomrule |
691 |
|
|
\end{tabular} |
692 |
gezelter |
3066 |
\label{tab:DielectricMap} |
693 |
chrisfen |
3064 |
\end{table} |
694 |
gezelter |
3066 |
|
695 |
chrisfen |
3064 |
The results of these calculations are displayed in table |
696 |
gezelter |
3066 |
\ref{tab:DielectricMap}. The dielectric constants for both models |
697 |
|
|
decrease linearly with increasing cutoff radii ($R_\textrm{c}$) and |
698 |
|
|
increase linearly with increasing damping ($\alpha$). Another point |
699 |
|
|
to note is that choosing $\alpha$ and $R_\textrm{c}$ identical to |
700 |
|
|
those used with the Ewald summation results in the same calculated |
701 |
|
|
dielectric constant. As an example, in the paper outlining the |
702 |
|
|
development of TIP5P-E, the real-space cutoff and Ewald coefficient |
703 |
|
|
were tethered to the system size, and for a 512 molecule system are |
704 |
|
|
approximately 12~\AA\ and 0.25~\AA$^{-1}$ respectively.\cite{Rick04} |
705 |
|
|
These parameters resulted in a dielectric constant of 92$\pm$14, while |
706 |
|
|
with SF these parameters give a dielectric constant of |
707 |
chrisfen |
3064 |
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where |
708 |
|
|
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and |
709 |
|
|
0.35~\AA$^{-1}$, and these parameters resulted in a dielectric |
710 |
gezelter |
3066 |
constant equal to 63$\pm$1.\cite{Horn04} Calculations using SF with |
711 |
|
|
these parameters and this water model give a dielectric constant of |
712 |
|
|
61$\pm$1. Since the dielectric constant is dependent on $\alpha$ and |
713 |
|
|
$R_\textrm{c}$ with the SF technique, it might be interesting to |
714 |
gezelter |
3068 |
investigate the dependence of the static dielectric constant on the |
715 |
|
|
choice of convergence parameters ($R_\textrm{c}$ and $\kappa$) |
716 |
|
|
utilized in most implementations of the Ewald sum. |
717 |
chrisfen |
3064 |
|
718 |
gezelter |
3068 |
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 |
gezelter |
3066 |
|
724 |
|
|
Although it is tempting to choose damping parameters equivalent to the |
725 |
gezelter |
3068 |
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 |
gezelter |
3066 |
|
746 |
gezelter |
3068 |
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 |
gezelter |
3066 |
|
753 |
gezelter |
3068 |
As a final note on optimal damping parameters, aside from a slight |
754 |
chrisfen |
3064 |
lowering of $\Delta H_\textrm{vap}$, using these $\alpha$ values does |
755 |
gezelter |
3068 |
not alter any of the other thermodynamic properties. |
756 |
chrisfen |
3064 |
|
757 |
gezelter |
3067 |
\subsection{Dynamic Properties of TIP5P-E}\label{sec:t5peDynamics} |
758 |
chrisfen |
3064 |
|
759 |
gezelter |
3068 |
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 |
chrisfen |
3064 |
|
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^\alpha(t) |
775 |
|
|
\cdot\hat{\mathbf{u}}_i^\alpha(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 |
gezelter |
3067 |
axis $\alpha$. 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 |
chrisfen |
3064 |
|
786 |
|
|
From the orientation autocorrelation functions, we can obtain time |
787 |
gezelter |
3067 |
constants for rotational relaxation. The relatively short time |
788 |
chrisfen |
3064 |
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 |
gezelter |
3067 |
\includegraphics[width=5.8in]{./figures/t5peDynamics.pdf} |
801 |
chrisfen |
3064 |
\caption{Diffusion constants ({\it upper}) and reorientational time |
802 |
gezelter |
3066 |
constants ({\it lower}) for TIP5P-E using the Ewald sum and SF |
803 |
chrisfen |
3064 |
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 |
gezelter |
3066 |
the Ewald sum are reproduced with the SF technique. The enhanced |
812 |
gezelter |
3067 |
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 |
gezelter |
3068 |
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 |
chrisfen |
3064 |
|
822 |
|
|
The $\tau_2^y$ results in the lower frame of figure |
823 |
gezelter |
3067 |
\ref{fig:t5peDynamics} show a much greater difference between the SF |
824 |
|
|
results and the Ewald results. At all temperatures shown, TIP5P-E |
825 |
chrisfen |
3064 |
relaxes faster than experiment with the Ewald sum while tracking |
826 |
gezelter |
3067 |
experiment fairly well when using the SF technique, independent of the |
827 |
gezelter |
3068 |
choice of damping constant. There are several possible reasons for |
828 |
gezelter |
3067 |
this deviation between techniques. The Ewald results were calculated |
829 |
gezelter |
3068 |
using shorter (10~ps) trajectories than the SF results (200~ps). |
830 |
|
|
Calculation of these SF values from a 10~ps trajectory (with |
831 |
|
|
subsequently lower accuracy) showed a 0.4~ps drop in $\tau_2^y$, |
832 |
gezelter |
3067 |
placing the result more in line with that obtained using the Ewald |
833 |
gezelter |
3068 |
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 |
chrisfen |
3064 |
|
850 |
gezelter |
3067 |
\subsection{Comparison of Reaction Field and Damped Electrostatics for |
851 |
|
|
SSD/RF} |
852 |
chrisfen |
3064 |
|
853 |
gezelter |
3068 |
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 |
chrisfen |
3064 |
|
864 |
gezelter |
3068 |
The properties shown in table \ref{tab:dampedSSDRF} show that damped |
865 |
|
|
electrostatics produces even better agreement with experiment than is |
866 |
|
|
obtained via reaction field. The average density shows a modest |
867 |
|
|
increase when using damped electrostatics in place of the reaction |
868 |
|
|
field. This comes about because we neglect the pressure effect due to |
869 |
|
|
the surroundings outside of the cutoff, instead relying on screening |
870 |
|
|
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. Using damped electrostatics, SSD/RF can be |
881 |
|
|
used effectively with mixed charge / dipolar systems, such as |
882 |
|
|
dissolved ions, dissolved organic molecules, or even proteins. |
883 |
chrisfen |
3064 |
|
884 |
|
|
\begin{table} |
885 |
gezelter |
3068 |
\caption{Properties of SSD/RF when using reaction field or damped |
886 |
|
|
electrostatics. Simulations were carried out at 298~K, 1~atm, and |
887 |
|
|
with 512 molecules.} |
888 |
chrisfen |
3064 |
\footnotesize |
889 |
|
|
\centering |
890 |
|
|
\begin{tabular}{ llccc } |
891 |
|
|
\toprule |
892 |
|
|
\toprule |
893 |
gezelter |
3068 |
& & Reaction Field (Ref. \citen{Fennell04}) & Damped Electrostatics & |
894 |
|
|
Experiment [Ref.] \\ |
895 |
chrisfen |
3064 |
& & $\epsilon = 80$ & $R_\textrm{c} = 12$\AA ; $\alpha = 0.2125$~\AA$^{-1}$ & \\ |
896 |
|
|
\midrule |
897 |
|
|
$\rho$ & (g cm$^{-3}$) & 0.997(1) & 1.004(4) & 0.997 [\citen{CRC80}]\\ |
898 |
|
|
$C_p$ & (cal mol$^{-1}$ K$^{-1}$) & 23.8(2) & 27(1) & 18.005 [\citen{Wagner02}] \\ |
899 |
|
|
$D$ & ($10^{-5}$ cm$^2$ s$^{-1}$) & 2.32(6) & 2.33(2) & 2.299 [\citen{Mills73}]\\ |
900 |
|
|
$n_C$ & & 4.4 & 4.2 & 4.7 [\citen{Hura00}]\\ |
901 |
|
|
$n_H$ & & 3.7 & 3.7 & 3.5 [\citen{Soper86}]\\ |
902 |
|
|
$\tau_1$ & (ps) & 7.2(4) & 5.86(8) & 5.7 [\citen{Eisenberg69}]\\ |
903 |
|
|
$\tau_2$ & (ps) & 3.2(2) & 2.45(7) & 2.3 [\citen{Krynicki66}]\\ |
904 |
|
|
\bottomrule |
905 |
|
|
\end{tabular} |
906 |
|
|
\label{tab:dampedSSDRF} |
907 |
|
|
\end{table} |
908 |
|
|
|
909 |
gezelter |
3067 |
\subsection{Predictions of Ice Polymorph Stability} |
910 |
chrisfen |
3064 |
|
911 |
gezelter |
3068 |
In an earlier paper, we performed a series of free energy calculations |
912 |
chrisfen |
3064 |
on several ice polymorphs which are stable or metastable at low |
913 |
|
|
pressures, one of which (Ice-$i$) we observed in spontaneous |
914 |
gezelter |
3068 |
crystallizations of an early version of the SSD/RF water |
915 |
|
|
model.\cite{Fennell05} In this study, a distinct dependence of the |
916 |
|
|
computed free energies on the cutoff radius and electrostatic |
917 |
|
|
summation method was observed. Since the SF technique can be used as |
918 |
|
|
a simple and efficient replacement for the Ewald summation, our final |
919 |
|
|
test of this method is to see if it is capable of addressing the |
920 |
|
|
spurious stability of the Ice-$i$ phases with the various common water |
921 |
|
|
models. To this end, we have performed thermodynamic integrations of |
922 |
|
|
all of the previously discussed ice polymorphs using the SF technique |
923 |
|
|
with a cutoff radius of 12~\AA\ and an $\alpha$ of 0.2125~\AA . These |
924 |
|
|
calculations were performed on TIP5P-E and TIP4P-Ew (variants of the |
925 |
|
|
TIP5P and TIP4P models optimized for the Ewald summation) as well as |
926 |
|
|
SPC/E and SSD/RF. |
927 |
chrisfen |
3064 |
|
928 |
|
|
\begin{table} |
929 |
|
|
\centering |
930 |
|
|
\caption{Helmholtz free energies of ice polymorphs at 1~atm and 200~K |
931 |
gezelter |
3066 |
using the damped SF electrostatic correction method with a |
932 |
chrisfen |
3064 |
variety of water models.} |
933 |
|
|
\begin{tabular}{ lccccc } |
934 |
|
|
\toprule |
935 |
|
|
\toprule |
936 |
|
|
Model & I$_\textrm{h}$ & I$_\textrm{c}$ & B & Ice-$i$ & Ice-$i^\prime$ \\ |
937 |
|
|
\cmidrule(lr){2-6} |
938 |
|
|
& \multicolumn{5}{c}{(kcal mol$^{-1}$)} \\ |
939 |
|
|
\midrule |
940 |
|
|
TIP5P-E & -11.98(4) & -11.96(4) & -11.87(3) & - & -11.95(3) \\ |
941 |
|
|
TIP4P-Ew & -13.11(3) & -13.09(3) & -12.97(3) & - & -12.98(3) \\ |
942 |
|
|
SPC/E & -12.99(3) & -13.00(3) & -13.03(3) & - & -12.99(3) \\ |
943 |
|
|
SSD/RF & -11.83(3) & -11.66(4) & -12.32(3) & -12.39(3) & - \\ |
944 |
|
|
\bottomrule |
945 |
|
|
\end{tabular} |
946 |
|
|
\label{tab:dampedFreeEnergy} |
947 |
|
|
\end{table} |
948 |
|
|
The results of these calculations in table \ref{tab:dampedFreeEnergy} |
949 |
gezelter |
3067 |
show similar behavior to the Ewald results in the previous |
950 |
|
|
study.\cite{Fennell05} The Helmholtz free energies of the ice |
951 |
|
|
polymorphs with SSD/RF order in the same fashion, with Ice-$i$ having |
952 |
|
|
the lowest free energy; however, the Ice-$i$ and ice B free energies |
953 |
|
|
are quite a bit closer (nearly isoenergetic). The SPC/E results show |
954 |
gezelter |
3068 |
the different polymorphs to be nearly isoenergetic. This is the same |
955 |
|
|
behavior observed using an Ewald correction.\cite{Fennell05} Ice B has |
956 |
|
|
the lowest Helmholtz free energy for SPC/E; however, all the polymorph |
957 |
|
|
results overlap within the error estimates. |
958 |
chrisfen |
3064 |
|
959 |
|
|
The most interesting results from these calculations come from the |
960 |
|
|
more expensive TIP4P-Ew and TIP5P-E results. Both of these models were |
961 |
|
|
optimized for use with an electrostatic correction and are |
962 |
gezelter |
3067 |
geometrically arranged to mimic water using drastically different |
963 |
gezelter |
3068 |
charge distributions. In TIP5P-E, the primary location for the |
964 |
gezelter |
3067 |
negative charge in the molecule is assigned to the lone-pairs of the |
965 |
|
|
oxygen, while TIP4P-Ew places the negative charge near the |
966 |
|
|
center-of-mass along the H-O-H bisector. There is some debate as to |
967 |
|
|
which is the proper choice for the negative charge location, and this |
968 |
|
|
has in part led to a six-site water model that balances both of these |
969 |
|
|
options.\cite{Vega05,Nada03} The limited results in table |
970 |
|
|
\ref{tab:dampedFreeEnergy} support the results of Vega {\it et al.}, |
971 |
gezelter |
3068 |
which indicate the TIP4P charge location geometry performs better |
972 |
|
|
under some circumstances.\cite{Vega05} With the TIP4P-Ew water model, |
973 |
|
|
the experimentally observed polymorph (ice I$_\textrm{h}$) is the |
974 |
|
|
preferred form with ice I$_\textrm{c}$ slightly higher in energy, |
975 |
|
|
though overlapping within error. Additionally, the spurious ice B and |
976 |
|
|
Ice-$i^\prime$ structures are destabilized relative to these |
977 |
|
|
polymorphs. TIP5P-E shows similar behavior to SPC/E, where there is no |
978 |
|
|
real free energy distinction between the various polymorphs. While ice |
979 |
|
|
B is close in free energy to the other polymorphs, these results fail |
980 |
|
|
to support the findings of other researchers indicating the preferred |
981 |
|
|
form of TIP5P at 1~atm is a structure similar to ice |
982 |
|
|
B.\cite{Yamada02,Vega05,Abascal05} It should be noted that we were |
983 |
|
|
looking at TIP5P-E rather than TIP5P, and the differences in the |
984 |
|
|
Lennard-Jones parameters could cause this discrepancy. Overall, these |
985 |
|
|
results indicate that TIP4P-Ew is a better mimic of the solid forms of |
986 |
|
|
water than some of the other models. |
987 |
chrisfen |
3064 |
|
988 |
|
|
\section{Conclusions} |
989 |
|
|
|
990 |
gezelter |
3067 |
This investigation of pairwise electrostatic summation techniques |
991 |
|
|
shows that there is a viable and computationally efficient alternative |
992 |
|
|
to the Ewald summation. The SF method (eq. \ref{eq:DSFPot}) has |
993 |
|
|
proven itself capable of reproducing structural, thermodynamic, and |
994 |
|
|
dynamic quantities that are nearly quantitative matches to results |
995 |
|
|
from far more expensive methods. Additionally, we have now extended |
996 |
|
|
the damping formalism to electrostatic multipoles, so the damped SF |
997 |
|
|
potential can be used in systems that contain mixtures of charges and |
998 |
|
|
point multipoles. |
999 |
|
|
|
1000 |
|
|
We have also provided a simple linear prescription for choosing |
1001 |
|
|
optimal damping parameters given a choice of cutoff radius. The |
1002 |
|
|
damping parameters were chosen to obtain a static dielectric constant |
1003 |
|
|
as close as possible to the experimental value, which should be useful |
1004 |
|
|
for simulating the electrostatic screening properties of liquid water |
1005 |
|
|
accurately. The linear formula for optimal damping was the same for |
1006 |
|
|
a complicated multipoint model as it was for a simple point-dipolar |
1007 |
|
|
model. |
1008 |
|
|
|
1009 |
|
|
As in all purely pairwise cutoff methods, the damped SF method is |
1010 |
|
|
expected to scale approximately {\it linearly} with system size, and |
1011 |
|
|
is easily parallelizable. This should result in substantial |
1012 |
|
|
reductions in the computational cost of performing large simulations. |
1013 |
|
|
With the proper use of pre-computation and spline interpolation, the |
1014 |
gezelter |
3068 |
damped SF method is essentially the same cost as a simple real-space |
1015 |
gezelter |
3067 |
cutoff. |
1016 |
|
|
|
1017 |
|
|
We are not suggesting that there is any flaw with the Ewald sum; in |
1018 |
|
|
fact, it is the standard by which the damped SF method has been |
1019 |
|
|
judged. However, these results provide further evidence that in the |
1020 |
|
|
typical simulations performed today, the Ewald summation may no longer |
1021 |
|
|
be required to obtain the level of accuracy most researchers have come |
1022 |
|
|
to expect. |
1023 |
|
|
|
1024 |
chrisfen |
3064 |
\section{Acknowledgments} |
1025 |
|
|
Support for this project was provided by the National Science |
1026 |
|
|
Foundation under grant CHE-0134881. Computation time was provided by |
1027 |
gezelter |
3067 |
the Notre Dame Center for Research Computing. The authors would like |
1028 |
|
|
to thank Steve Corcelli and Ed Maginn for helpful discussions and |
1029 |
|
|
comments. |
1030 |
chrisfen |
3064 |
|
1031 |
|
|
\newpage |
1032 |
|
|
|
1033 |
|
|
\bibliographystyle{achemso} |
1034 |
|
|
\bibliography{multipoleSFPaper} |
1035 |
|
|
|
1036 |
|
|
|
1037 |
|
|
\end{document} |