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