| 1 |
|
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
| 2 |
< |
\documentclass[12pt]{article} |
| 2 |
> |
\documentclass[11pt]{article} |
| 3 |
|
\usepackage{tabls} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
|
\usepackage[tbtags]{amsmath} |
| 20 |
|
\brokenpenalty=10000 |
| 21 |
|
\renewcommand{\baselinestretch}{1.2} |
| 22 |
|
\renewcommand\citemid{\ } % no comma in optional reference note |
| 23 |
+ |
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions |
| 24 |
+ |
\let\Caption\caption |
| 25 |
+ |
\renewcommand\caption[1]{% |
| 26 |
+ |
\Caption[#1]{}% |
| 27 |
+ |
} |
| 28 |
+ |
\setlength{\abovecaptionskip}{1.2in} |
| 29 |
+ |
\setlength{\belowcaptionskip}{1.2in} |
| 30 |
|
|
| 31 |
|
\begin{document} |
| 32 |
|
|
| 33 |
< |
\title{Pairwise electrostatic corrections for monopolar and multipolar systems} |
| 33 |
> |
\title{Pairwise Alternatives to the Ewald Sum: Applications |
| 34 |
> |
and Extension to Point Multipoles} |
| 35 |
|
|
| 36 |
|
\author{Christopher J. Fennell and J. Daniel Gezelter \\ |
| 37 |
|
Department of Chemistry and Biochemistry\\ University of Notre Dame\\ |
| 43 |
|
%\doublespacing |
| 44 |
|
|
| 45 |
|
\begin{abstract} |
| 46 |
+ |
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 |
|
\end{abstract} |
| 57 |
|
|
| 58 |
|
%\narrowtext |
| 63 |
|
|
| 64 |
|
\section{Introduction} |
| 65 |
|
|
| 66 |
< |
Over the past several years, there has been increased interest in |
| 67 |
< |
pairwise methods for correcting electrostatic interaction in computer |
| 66 |
> |
Over the past several years, there has been increasing interest in |
| 67 |
> |
pairwise methods for correcting electrostatic interactions in computer |
| 68 |
|
simulations of condensed molecular |
| 69 |
< |
systems.\cite{Wolf99,Zahn02,Kast03,Fennell06} |
| 69 |
> |
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 |
|
|
| 76 |
< |
\section{Shifted Force Pairwise Electrostatic Correction} |
| 77 |
< |
|
| 78 |
< |
In our previous look at pairwise electrostatic correction methods, we |
| 79 |
< |
described the undamped and damped shifted potential {\sc sp} and |
| 80 |
< |
shifted force {\sc sf} techniques.\cite{Fennell06} These techniques |
| 81 |
< |
were developed from the observations and efforts of Wolf {\it et al.} |
| 82 |
< |
and their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99} |
| 83 |
< |
The potential for the damped form of the {\sc sp} method is given by |
| 76 |
> |
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 |
|
\begin{equation} |
| 62 |
– |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r} |
| 63 |
– |
- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) |
| 64 |
– |
\quad r\leqslant R_\textrm{c}, |
| 65 |
– |
\label{eq:DSPPot} |
| 66 |
– |
\end{equation} |
| 67 |
– |
while the damped form of the {\sc sf} method is given by |
| 68 |
– |
\begin{equation} |
| 85 |
|
\begin{split} |
| 86 |
|
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& |
| 87 |
|
\frac{\mathrm{erfc}\left(\alpha r\right)}{r} |
| 94 |
|
\label{eq:DSFPot} |
| 95 |
|
\end{split} |
| 96 |
|
\end{equation} |
| 97 |
< |
Overall, the damped {\sc sf} method is considered an improvement over |
| 98 |
< |
the {\sc sp} method because there is no discontinuity in the forces as |
| 83 |
< |
particles move out of the cutoff radius ($R_\textrm{c}$). This is |
| 84 |
< |
accomplished by shifting the forces (and potential) to zero at |
| 85 |
< |
$R_\textrm{c}$. This is analogous to neutralizing the charge and the |
| 86 |
< |
force effect of the charges within $R_\textrm{c}$. |
| 97 |
> |
(In this potential and in all electrostatic quantities that follow, |
| 98 |
> |
the standard $4 \pi \epsilon_{0}$ has been omitted for clarity.) |
| 99 |
|
|
| 100 |
< |
As we mentioned in the previous study, to complete the charge |
| 101 |
< |
neutralization procedure, a self-neutralization term needs to be |
| 102 |
< |
included in the potential. This term exists outside the pair-loop, and |
| 103 |
< |
can be thought of as a charge opposite in sign and equal in magnitude |
| 104 |
< |
to the central particle, spread out over the entire cutoff |
| 105 |
< |
sphere. This term is calculated via a single loop over all charges in |
| 106 |
< |
the system. For the undamped case, the self-term can be written as |
| 100 |
> |
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 |
|
\begin{equation} |
| 116 |
< |
V_\textrm{self} = \sum_{i=1}^N\frac{q_iq_i}{2R_\textrm{c}}, |
| 116 |
> |
V_\textrm{self} = \frac{1}{2 R_\textrm{c}} \sum_{i=1}^N q_i^{2}, |
| 117 |
|
\label{eq:selfTerm} |
| 118 |
|
\end{equation} |
| 119 |
|
while for the damped case it can be written as |
| 120 |
|
\begin{equation} |
| 121 |
< |
V_\textrm{Dself} = \sum_{i=1}^Nq_iq_i\left(\frac{\alpha}{\sqrt{\pi}} |
| 122 |
< |
+ \frac{\textrm{erfc}(\alpha R_\textrm{c})}{2R_\textrm{c}}\right). |
| 121 |
> |
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 |
|
\label{eq:dampSelfTerm} |
| 125 |
|
\end{equation} |
| 126 |
|
The first term within the parentheses of equation |
| 127 |
< |
(\ref{eq:dampSelfTerm}) is identical to the self-term in the Ewald |
| 127 |
> |
(\ref{eq:dampSelfTerm}) is identical to the self term in the Ewald |
| 128 |
|
summation, and comes from the utilization of the complimentary error |
| 129 |
|
function for electrostatic damping.\cite{deLeeuw80,Wolf99} |
| 130 |
|
|
| 131 |
< |
\section{Multipolar Electrostatic Damping}\label{sec:dampingMultipoles} |
| 131 |
> |
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 |
|
|
| 141 |
< |
As stated above, the {\sc sf} method operates by neutralizing the |
| 142 |
< |
charge pair force interactions (and potential) within the cutoff |
| 114 |
< |
sphere with shifting and by damping the electrostatic |
| 115 |
< |
interactions. Now we would like to consider an extension of these |
| 116 |
< |
techniques to include point multipole interactions. How will the |
| 117 |
< |
shifting and damping need to be modified in order to accommodate point |
| 118 |
< |
multipoles? |
| 141 |
> |
\section{Electrostatic Damping for Point |
| 142 |
> |
Multipoles}\label{sec:dampingMultipoles} |
| 143 |
|
|
| 144 |
< |
Of the two component techniques, the easiest to adapt is |
| 145 |
< |
shifting. Shifting is employed to neutralize the cutoff sphere; |
| 146 |
< |
however, in a system composed purely of point multipoles, the cutoff |
| 123 |
< |
sphere is already neutralized. This means that shifting is not |
| 124 |
< |
necessary between point multipoles. In a mixed system of monopoles and |
| 125 |
< |
multipoles, the undamped {\sc sf} potential needs only to shift the |
| 126 |
< |
force terms of the monopole and smoothly truncate the multipole |
| 127 |
< |
interactions with a switching function. The switching function is |
| 128 |
< |
required in order to conserve energy, because a discontinuity will |
| 129 |
< |
exist in both the potential and forces at $R_\textrm{c}$ in the |
| 130 |
< |
absence of shifting terms. |
| 144 |
> |
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 |
|
|
| 148 |
< |
If we want to dampen the {\sc sf} potential, then we need to |
| 149 |
< |
incorporate the complimentary error function term into the multipole |
| 150 |
< |
potentials. The most direct way to do this is by replacing $r^{-1}$ |
| 151 |
< |
with erfc$(\alpha r)\cdot r^{-1}$ in the multipole expansion. In the |
| 152 |
< |
multipole expansion, rather than considering only the interactions |
| 153 |
< |
between single point charges, the electrostatic interaction is |
| 154 |
< |
reformulated such that it describes the interaction between charge |
| 155 |
< |
distributions about central sites of the respective sets of |
| 156 |
< |
charges.\cite{Hirschfelder67} This procedure is what leads to the |
| 157 |
< |
familiar charge-dipole, |
| 148 |
> |
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 |
|
\begin{equation} |
| 175 |
< |
V_\textrm{cd} = -q_i\frac{\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}}{r^3_{ij}} |
| 176 |
< |
= -q_i\mu_j\frac{\cos\theta}{r^2_{ij}}, |
| 145 |
< |
\label{eq:chargeDipole} |
| 175 |
> |
T = \frac{1}{r_{ij}}, |
| 176 |
> |
\label{eq:tensorRank1} |
| 177 |
|
\end{equation} |
| 147 |
– |
and dipole-dipole, |
| 178 |
|
\begin{equation} |
| 179 |
< |
V_\textrm{dd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}) |
| 180 |
< |
(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}} - |
| 151 |
< |
\frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}}, |
| 152 |
< |
\label{eq:dipoleDipole} |
| 179 |
> |
T_\alpha = \nabla_\alpha \frac{1}{r_{ij}}, |
| 180 |
> |
\label{eq:tensorRank2} |
| 181 |
|
\end{equation} |
| 182 |
< |
interaction potentials. |
| 182 |
> |
\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 |
|
|
| 192 |
< |
Using the charge-dipole interaction as an example, if we insert |
| 193 |
< |
erfc$(\alpha r)\cdot r^{-1}$ in place of $r^{-1}$, a damped |
| 194 |
< |
charge-dipole results, |
| 192 |
> |
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 |
|
\begin{equation} |
| 197 |
|
V_\textrm{Dcd} = -q_i\mu_j\frac{\cos\theta}{r^2_{ij}} c_1(r_{ij}), |
| 198 |
|
\label{eq:dChargeDipole} |
| 203 |
|
+ \textrm{erfc}(\alpha r_{ij}). |
| 204 |
|
\label{eq:c1Func} |
| 205 |
|
\end{equation} |
| 169 |
– |
Thus, $c_1(r_{ij})$ is the resulting damping term that modifies the |
| 170 |
– |
standard charge-dipole potential (Eq. (\ref{eq:chargeDipole})). Note |
| 171 |
– |
that this damping term is dependent upon distance and not upon |
| 172 |
– |
orientation, and that it is acting on what was originally an |
| 173 |
– |
$r^{-3}$ function. By writing the damped form in this manner, we |
| 174 |
– |
can collect the damping into one function and apply it to the original |
| 175 |
– |
potential when damping is desired. This works well for potentials that |
| 176 |
– |
have only one $r^{-n}$ term (where $n$ is an odd positive integer); |
| 177 |
– |
but in the case of the dipole-dipole potential, there is one part |
| 178 |
– |
dependent on $r^{-3}$ and another dependent on $r^{-5}$. In order to |
| 179 |
– |
properly damping this potential, each of these parts is dampened with |
| 180 |
– |
separate damping functions. We can determine the necessary damping |
| 181 |
– |
functions by continuing with the multipole expansion; however, it |
| 182 |
– |
quickly becomes more complex with ``two-center'' systems, like the |
| 183 |
– |
dipole-dipole potential, and is typically approached with a spherical |
| 184 |
– |
harmonic formalism.\cite{Hirschfelder67} A simpler method for |
| 185 |
– |
determining these functions arises from adopting the tensor formalism |
| 186 |
– |
for expressing the electrostatic interactions.\cite{Stone02} |
| 206 |
|
|
| 207 |
< |
The tensor formalism for electrostatic interactions involves obtaining |
| 189 |
< |
the multipole interactions from successive gradients of the monopole |
| 190 |
< |
potential. Thus, tensors of rank one through three are |
| 191 |
< |
\begin{equation} |
| 192 |
< |
T = \frac{1}{4\pi\epsilon_0r_{ij}}, |
| 193 |
< |
\label{eq:tensorRank1} |
| 194 |
< |
\end{equation} |
| 207 |
> |
Equation \ref{eq:tensorRank3} generates a damped dipole-dipole potential, |
| 208 |
|
\begin{equation} |
| 209 |
< |
T_\alpha = \frac{1}{4\pi\epsilon_0}\nabla_\alpha \frac{1}{r_{ij}}, |
| 210 |
< |
\label{eq:tensorRank2} |
| 209 |
> |
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 |
|
\end{equation} |
| 216 |
+ |
where |
| 217 |
|
\begin{equation} |
| 200 |
– |
T_{\alpha\beta} = \frac{1}{4\pi\epsilon_0} |
| 201 |
– |
\nabla_\alpha\nabla_\beta \frac{1}{r_{ij}}, |
| 202 |
– |
\label{eq:tensorRank3} |
| 203 |
– |
\end{equation} |
| 204 |
– |
where the form of the first tensor gives the monopole-monopole |
| 205 |
– |
potential, the second gives the monopole-dipole potential, and the |
| 206 |
– |
third gives the monopole-quadrupole and dipole-dipole |
| 207 |
– |
potentials.\cite{Stone02} Since the force is $-\nabla V$, the forces |
| 208 |
– |
for each potential come from the next higher tensor. |
| 209 |
– |
|
| 210 |
– |
To obtain the damped electrostatic forms, we replace $r^{-1}$ with |
| 211 |
– |
erfc$(\alpha r)\cdot r^{-1}$. Equation \ref{eq:tensorRank2} generates |
| 212 |
– |
$c_1(r_{ij})$, just like the multipole expansion, while equation |
| 213 |
– |
\ref{eq:tensorRank3} generates $c_2(r_{ij})$, where |
| 214 |
– |
\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}). |
| 240 |
|
\label{eq:doubleFactorial} |
| 241 |
|
\end{equation} |
| 242 |
|
and $c_0(r_{ij})$ is erfc$(\alpha r_{ij})$. This generating function |
| 243 |
< |
is similar in form to those obtained by researchers for the |
| 244 |
< |
application of the Ewald sum to |
| 243 |
> |
is similar in form to those obtained by Smith and Aguado and Madden |
| 244 |
> |
for the application of the Ewald sum to |
| 245 |
|
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 |
< |
$r^{-3}$. In the damped dipole-dipole potential, |
| 247 |
< |
\begin{equation} |
| 248 |
< |
V_\textrm{Ddd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}) |
| 249 |
< |
(\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}} |
| 250 |
< |
c_2(r_{ij}) - |
| 251 |
< |
\frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}} |
| 252 |
< |
c_1(r_{ij}), |
| 253 |
< |
\label{eq:dampDipoleDipole} |
| 254 |
< |
\end{equation} |
| 249 |
> |
$r^{-3}$. |
| 250 |
|
$c_2(r_{ij})$ and $c_1(r_{ij})$ dampen these two parts |
| 251 |
< |
respectively. The forces for the damped dipole-dipole interaction, |
| 251 |
> |
respectively. The forces for the damped dipole-dipole interaction, are |
| 252 |
> |
obtained from the next higher tensor, $T_{\alpha \beta \gamma}$, |
| 253 |
|
\begin{equation} |
| 254 |
|
\begin{split} |
| 255 |
|
F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}) |
| 262 |
|
\end{split} |
| 263 |
|
\label{eq:dampDipoleDipoleForces} |
| 264 |
|
\end{equation} |
| 265 |
< |
rely on higher order damping functions because we perform another |
| 266 |
< |
gradient operation. In this manner, we can dampen higher order |
| 267 |
< |
multipolar interactions along with the monopole interactions, allowing |
| 268 |
< |
us to include multipoles in simulations involving damped electrostatic |
| 273 |
< |
interactions. |
| 265 |
> |
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 |
|
|
| 270 |
< |
\section{Applications of Pairwise Electrostatic Corrections: Water} |
| 270 |
> |
\section{Applications of DSF Electrostatics} |
| 271 |
|
|
| 272 |
< |
In our earlier work, we performed a statistical comparison of both |
| 273 |
< |
energetics and dynamics of the {\sc sf} and {\sc sp} methods, showing |
| 274 |
< |
that these pairwise correction techniques are nearly equivalent to the |
| 275 |
< |
Ewald summation in the simulation of general condensed |
| 276 |
< |
systems.\cite{Fennell06} It would be beneficial to have some specific |
| 277 |
< |
examples to better illustrate the similarities and differences of the |
| 278 |
< |
pairwise (specifically {\sc sf}) and Ewald techniques. To address |
| 279 |
< |
this, we decided to perform a detailed analysis involving water |
| 280 |
< |
simulations. Water is a good test in that there is a great deal of |
| 281 |
< |
detailed experimental information as well as an extensive library of |
| 282 |
< |
results from a variety of theoretical models. In choosing a model for |
| 288 |
< |
comparing these correction techniques, it is important to consider a |
| 289 |
< |
model that has been optimized for use with corrected electrostatics so |
| 290 |
< |
that the calculated properties are in better agreement with |
| 291 |
< |
experiment.\cite{vanderSpoel98,Horn04,Rick04} For this reason, we |
| 292 |
< |
chose the TIP5P-E water model.\cite{Rick04} |
| 272 |
> |
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 |
|
|
| 284 |
< |
\subsection{TIP5P-E} |
| 284 |
> |
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 |
|
|
| 293 |
+ |
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 |
|
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 |
< |
maximum anomaly present in liquid water near 4$^\circ$C. As with many |
| 316 |
< |
previous point charge water models (such as ST2, TIP3P, TIP4P, SPC, |
| 317 |
< |
and SPC/E), TIP5P was parametrized using a simple cutoff with no |
| 318 |
< |
long-range electrostatic |
| 315 |
> |
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 |
|
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 |
< |
conditions. When using any form of long-range electrostatic |
| 325 |
< |
correction, it has become common practice to develop or utilize a |
| 326 |
< |
reparametrized water model that corrects for this |
| 311 |
< |
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows |
| 312 |
< |
this practice and was optimized for use with the Ewald |
| 313 |
< |
summation.\cite{Rick04} In his publication, Rick preserved the |
| 314 |
< |
geometry and point charge magnitudes in TIP5P and focused on altering |
| 315 |
< |
the Lennard-Jones parameters to correct the density at 298~K. With the |
| 324 |
> |
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 |
|
density corrected, he compared common water properties for TIP5P-E |
| 328 |
|
using the Ewald sum with TIP5P using a 9~\AA\ cutoff. |
| 329 |
|
|
| 330 |
< |
In the following sections, we compared these same water properties |
| 331 |
< |
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 332 |
< |
{\sc sf} technique. In the above evaluation of the pairwise |
| 333 |
< |
techniques, we observed some flexibility in the choice of parameters. |
| 334 |
< |
Because of this, the following comparisons include the {\sc sf} |
| 324 |
< |
technique with a 12~\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and |
| 325 |
< |
0.2~\AA$^{-1}$, as well as a 9~\AA\ cutoff with an $\alpha$ = |
| 326 |
< |
0.2~\AA$^{-1}$. |
| 330 |
> |
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 |
|
|
| 336 |
< |
\subsubsection{Density}\label{sec:t5peDensity} |
| 336 |
> |
\subsection{The Density Maximum of TIP5P-E}\label{sec:t5peDensity} |
| 337 |
|
|
| 330 |
– |
As stated previously, the property that prompted the development of |
| 331 |
– |
TIP5P-E was the density at 1 atm. The density depends upon the |
| 332 |
– |
internal pressure of the system in the $NPT$ ensemble, and the |
| 333 |
– |
calculation of the pressure includes a components from both the |
| 334 |
– |
kinetic energy and the virial. More specifically, the instantaneous |
| 335 |
– |
molecular pressure ($p(t)$) is given by |
| 336 |
– |
\begin{equation} |
| 337 |
– |
p(t) = \frac{1}{\textrm{d}V}\sum_\mu |
| 338 |
– |
\left[\frac{\mathbf{P}_{\mu}^2}{M_{\mu}} |
| 339 |
– |
+ \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right], |
| 340 |
– |
\label{eq:MolecularPressure} |
| 341 |
– |
\end{equation} |
| 342 |
– |
where d is the dimensionality of the system, $V$ is the volume, |
| 343 |
– |
$\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$ |
| 344 |
– |
is the position of the center of mass ($M_\mu$) of molecule $\mu$, and |
| 345 |
– |
$\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule |
| 346 |
– |
$\mu$.\cite{Melchionna93} The virial term (the right term in the |
| 347 |
– |
brackets of equation \ref{eq:MolecularPressure}) is directly dependent |
| 348 |
– |
on the interatomic forces. Since the {\sc sp} method does not modify |
| 349 |
– |
the forces, the pressure using {\sc sp} will be identical to that |
| 350 |
– |
obtained without an electrostatic correction. The {\sc sf} method |
| 351 |
– |
does alter the virial component and, by way of the modified pressures, |
| 352 |
– |
should provide densities more in line with those obtained using the |
| 353 |
– |
Ewald summation. |
| 354 |
– |
|
| 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 |
| 350 |
|
\begin{figure} |
| 351 |
|
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf} |
| 352 |
|
\caption{Density versus temperature for the TIP5P-E water model when |
| 353 |
< |
using the Ewald summation (Ref. \citen{Rick04}) and the {\sc sf} method |
| 354 |
< |
with various parameters. The pressure term from the image-charge shell |
| 355 |
< |
is larger than that provided by the reciprocal-space portion of the |
| 356 |
< |
Ewald summation, leading to slightly lower densities. This effect is |
| 357 |
< |
more visible with the 9~\AA\ cutoff, where the image charges exert a |
| 358 |
< |
greater force on the central particle. The error bars for the {\sc sf} |
| 359 |
< |
methods show the average one-sigma uncertainty of the density |
| 360 |
< |
measurement, and this uncertainty is the same for all the {\sc sf} |
| 361 |
< |
curves.} |
| 353 |
> |
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 |
|
\label{fig:t5peDensities} |
| 363 |
|
\end{figure} |
| 364 |
|
Figure \ref{fig:t5peDensities} shows the densities calculated for |
| 365 |
< |
TIP5P-E using differing electrostatic corrections overlaid on the |
| 366 |
< |
experimental values.\cite{CRC80} The densities when using the {\sc sf} |
| 367 |
< |
technique are close to, though typically lower than, those calculated |
| 365 |
> |
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 |
|
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 |
< |
Ewald summation. Bringing the image charges closer to the central |
| 371 |
> |
Ewald summation. Bringing the image charges closer to the central |
| 372 |
|
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 |
| 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 |
< |
the damped {\sc sf} sums showing slightly higher densities; however, |
| 385 |
< |
it is apparent that the choice of cutoff radius plays a much more |
| 386 |
< |
important role in the resulting densities. |
| 384 |
> |
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 |
|
|
| 388 |
< |
As a final note, all of the above density calculations were performed |
| 389 |
< |
with systems of 512 water molecules. Rick observed a system size |
| 390 |
< |
dependence of the computed densities when using the Ewald summation, |
| 391 |
< |
most likely due to his tying of the convergence parameter to the box |
| 388 |
> |
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 |
|
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 |
< |
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also |
| 396 |
< |
lower the densities. Moving to larger systems, as long as the |
| 397 |
< |
R$_\textrm{c}$ remains at a fixed value, we would expect the densities |
| 398 |
< |
to remain constant. |
| 395 |
> |
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 |
|
|
| 400 |
< |
\subsubsection{Liquid Structure}\label{sec:t5peLiqStructure} |
| 400 |
> |
\subsection{Liquid Structure of TIP5P-E}\label{sec:t5peLiqStructure} |
| 401 |
|
|
| 419 |
– |
A common function considered when developing and comparing water |
| 420 |
– |
models is the oxygen-oxygen radial distribution function |
| 421 |
– |
($g_\textrm{OO}(r)$). The $g_\textrm{OO}(r)$ is the probability of |
| 422 |
– |
finding a pair of oxygen atoms some distance ($r$) apart relative to a |
| 423 |
– |
random distribution at the same density.\cite{Allen87} It is |
| 424 |
– |
calculated via |
| 425 |
– |
\begin{equation} |
| 426 |
– |
g_\textrm{OO}(r) = \frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i} |
| 427 |
– |
\delta(\mathbf{r}-\mathbf{r}_{ij})\right\rangle, |
| 428 |
– |
\label{eq:GOOofR} |
| 429 |
– |
\end{equation} |
| 430 |
– |
where the double sum is over all $i$ and $j$ pairs of $N$ oxygen |
| 431 |
– |
atoms. The $g_\textrm{OO}(r)$ can be directly compared to X-ray or |
| 432 |
– |
neutron scattering experiments through the oxygen-oxygen structure |
| 433 |
– |
factor ($S_\textrm{OO}(k)$) by the following relationship: |
| 434 |
– |
\begin{equation} |
| 435 |
– |
S_\textrm{OO}(k) = 1 + 4\pi\rho |
| 436 |
– |
\int_0^\infty r^2\frac{\sin kr}{kr}g_\textrm{OO}(r)\textrm{d}r. |
| 437 |
– |
\label{eq:SOOofK} |
| 438 |
– |
\end{equation} |
| 439 |
– |
Thus, $S_\textrm{OO}(k)$ is related to the Fourier-Laplace transform |
| 440 |
– |
of $g_\textrm{OO}(r)$. |
| 441 |
– |
|
| 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 |
< |
check whether the choice of using the Ewald summation or the {\sc sf} |
| 407 |
> |
check whether the choice of using the Ewald summation or the SF |
| 408 |
|
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. |
| 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 |
< |
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially |
| 417 |
> |
using the SF technique, the $g_\textrm{OO}(r)$s are essentially |
| 418 |
|
identical.} |
| 419 |
|
\label{fig:t5peGofRs} |
| 420 |
|
\end{figure} |
| 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 |
< |
are indistinguishable. These results indicate that the |
| 426 |
> |
are indistinguishable. These results do indicate that |
| 427 |
|
$g_\textrm{OO}(r)$ is insensitive to the choice of electrostatic |
| 428 |
|
correction. |
| 429 |
|
|
| 430 |
< |
\subsubsection{Thermodynamic Properties}\label{sec:t5peThermo} |
| 430 |
> |
\subsection{Thermodynamic Properties of TIP5P-E}\label{sec:t5peThermo} |
| 431 |
|
|
| 432 |
< |
In addition to the density, there are a variety of thermodynamic |
| 433 |
< |
quantities that can be calculated for water and compared directly to |
| 434 |
< |
experimental values. Some of these additional quantities include the |
| 435 |
< |
latent heat of vaporization ($\Delta H_\textrm{vap}$), the constant |
| 436 |
< |
pressure heat capacity ($C_p$), the isothermal compressibility |
| 437 |
< |
($\kappa_T$), the thermal expansivity ($\alpha_p$), and the static |
| 438 |
< |
dielectric constant ($\epsilon$). All of these properties were |
| 439 |
< |
calculated for TIP5P-E with the Ewald summation, so they provide a |
| 440 |
< |
good set for comparisons involving the {\sc sf} technique. |
| 432 |
> |
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 |
|
|
| 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 |
| 445 |
|
determined via |
| 446 |
|
\begin{equation} |
| 447 |
|
\begin{split} |
| 448 |
< |
\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, |
| 448 |
> |
\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 |
|
\end{split} |
| 453 |
|
\label{eq:DeltaHVap} |
| 454 |
|
\end{equation} |
| 455 |
< |
where $E$ is the total energy, $U$ is the potential energy, $p$ is the |
| 455 |
> |
where $E$ is the total energy, $U$ is the potential energy, $P$ is the |
| 456 |
|
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 |
< |
the $U_\textrm{liq.}$ term. Additionally, the $pV$ term for the gas is |
| 461 |
> |
the $U_\textrm{liq}$ term. Additionally, the $pV$ term for the gas is |
| 462 |
|
several orders of magnitude larger than that of the liquid, so we can |
| 463 |
|
neglect the liquid $pV$ term. |
| 464 |
|
|
| 468 |
|
enthalpy in constant pressure simulations via |
| 469 |
|
\begin{equation} |
| 470 |
|
\begin{split} |
| 471 |
< |
C_p = \left(\frac{\partial H}{\partial T}\right)_{N,p} |
| 471 |
> |
C_p = \left(\frac{\partial H}{\partial T}\right)_{N,P} |
| 472 |
|
= \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2}, |
| 473 |
|
\end{split} |
| 474 |
|
\label{eq:Cp} |
| 477 |
|
\begin{equation} |
| 478 |
|
\begin{split} |
| 479 |
|
\kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T} |
| 480 |
< |
= \frac{(\langle V^2\rangle_{N,P,T} - \langle V\rangle^2_{N,P,T})} |
| 481 |
< |
{k_BT\langle V\rangle_{N,P,T}}, |
| 480 |
> |
= \frac{(\langle V^2\rangle_{NPT} - \langle V\rangle^2_{NPT})} |
| 481 |
> |
{k_BT\langle V\rangle_{NPT}}, |
| 482 |
|
\end{split} |
| 483 |
|
\label{eq:kappa} |
| 484 |
|
\end{equation} |
| 486 |
|
\begin{equation} |
| 487 |
|
\begin{split} |
| 488 |
|
\alpha_p = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P} |
| 489 |
< |
= \frac{(\langle VH\rangle_{N,P,T} |
| 490 |
< |
- \langle V\rangle_{N,P,T}\langle H\rangle_{N,P,T})} |
| 491 |
< |
{k_BT^2\langle V\rangle_{N,P,T}}. |
| 489 |
> |
= \frac{(\langle VH\rangle_{NPT} |
| 490 |
> |
- \langle V\rangle_{NPT}\langle H\rangle_{NPT})} |
| 491 |
> |
{k_BT^2\langle V\rangle_{NPT}}. |
| 492 |
|
\end{split} |
| 493 |
|
\label{eq:alpha} |
| 494 |
|
\end{equation} |
| 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 |
< |
equal to zero. Since the {\sc sf} method also lacks this |
| 529 |
> |
equal to zero. Since the SF method also lacks this |
| 530 |
|
dielectric/surface term, equation (\ref{eq:staticDielectric}) is still |
| 531 |
|
valid for determining static dielectric constants. |
| 532 |
|
|
| 541 |
|
\centering |
| 542 |
|
\includegraphics[width=4.5in]{./figures/t5peThermo.pdf} |
| 543 |
|
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation |
| 544 |
< |
and the {\sc sf} techniques along with the experimental values. Units |
| 544 |
> |
and the SF techniques along with the experimental values. Units |
| 545 |
|
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 |
| 553 |
|
\label{fig:t5peThermo} |
| 554 |
|
\end{figure} |
| 555 |
|
|
| 556 |
< |
As observed for the density in section \ref{sec:t5peDensity}, the |
| 557 |
< |
property trends with temperature seen when using the Ewald summation |
| 558 |
< |
are reproduced with the {\sc sf} technique. One noticeable difference |
| 559 |
< |
between the properties calculated using the two methods are the lower |
| 560 |
< |
$\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be |
| 561 |
< |
expected due to the direct weakening of the electrostatic interaction |
| 562 |
< |
through forced neutralization. This results in an increase of the |
| 563 |
< |
intermolecular potential producing lower values from equation |
| 564 |
< |
(\ref{eq:DeltaHVap}). The slopes of these values with temperature are |
| 565 |
< |
similar to that seen using the Ewald summation; however, they are both |
| 566 |
< |
steeper than the experimental trend, indirectly resulting in the |
| 607 |
< |
inflated $C_p$ values at all temperatures. |
| 556 |
> |
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 |
|
|
| 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 |
< |
and $C_p$ results discussed in the previous paragraph, the deviations |
| 571 |
< |
between experiment and simulation in this region are not the fault of |
| 572 |
< |
the electrostatic summation methods but are due to the geometry and |
| 573 |
< |
parameters of the TIP5P class of water models. Like most rigid, |
| 574 |
< |
non-polarizable, point-charge water models, the density decreases with |
| 575 |
< |
temperature at a much faster rate than experiment (see figure |
| 576 |
< |
\ref{fig:t5peDensities}). This reduced density leads to the inflated |
| 577 |
< |
compressibility and expansivity values at higher temperatures seen |
| 578 |
< |
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability |
| 579 |
< |
and many-body effects are required in order for water models to |
| 580 |
< |
overcome differences between simulation-based and experimentally |
| 622 |
< |
determined densities at these higher |
| 570 |
> |
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 |
|
temperatures.\cite{Laasonen93,Donchev06} |
| 582 |
|
|
| 583 |
|
At temperatures below the freezing point for experimental water, the |
| 584 |
< |
differences between {\sc sf} and the Ewald summation results are more |
| 584 |
> |
differences between SF and the Ewald summation results are more |
| 585 |
|
apparent. The larger $C_p$ and lower $\alpha_p$ values in this region |
| 586 |
|
indicate a more pronounced transition in the supercooled regime, |
| 587 |
< |
particularly in the case of {\sc sf} without damping. This points to |
| 630 |
< |
the onset of a more frustrated or glassy behavior for TIP5P-E at |
| 631 |
< |
temperatures below 250~K in the {\sc sf} simulations, indicating that |
| 632 |
< |
disorder in the reciprocal-space term of the Ewald summation might act |
| 633 |
< |
to loosen up the local structure more than the image-charges in {\sc |
| 634 |
< |
sf}. The damped {\sc sf} actually makes a better comparison with |
| 635 |
< |
experiment in this region, particularly for the $\alpha_p$ values. The |
| 636 |
< |
local interactions in the undamped {\sc sf} technique appear to be too |
| 637 |
< |
strong since the property change is much more dramatic than the damped |
| 638 |
< |
forms, while the Ewald summation appears to weight the |
| 639 |
< |
reciprocal-space interactions at the expense the local interactions, |
| 640 |
< |
disagreeing with the experimental results. This observation is |
| 641 |
< |
explored further in section \ref{sec:t5peDynamics}. |
| 587 |
> |
particularly in the case of SF without damping. |
| 588 |
|
|
| 589 |
+ |
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 |
|
The final thermodynamic property displayed in figure |
| 599 |
|
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy |
| 600 |
< |
between the Ewald summation and the {\sc sf} technique (and experiment |
| 601 |
< |
for that matter). It is known that the dielectric constant is |
| 602 |
< |
dependent upon and quite sensitive to the imposed boundary |
| 603 |
< |
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the |
| 604 |
< |
converged $\epsilon$ values accumulated for the {\sc sf} |
| 650 |
< |
simulations. Lack of a damping function results in dielectric |
| 600 |
> |
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 |
|
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 |
< |
``Ewald coefficient'' value also has a significant effect on the |
| 609 |
< |
calculated value when using the Ewald summation. In the simulations of |
| 610 |
< |
TIP5P-E with the Ewald sum, this screening parameter was tethered to |
| 611 |
< |
the simulation box size (as was the $R_\textrm{c}$).\cite{Rick04} In |
| 612 |
< |
general, systems with larger screening parameters reported larger |
| 613 |
< |
dielectric constant values, the same behavior we see here with {\sc |
| 614 |
< |
sf}; however, the choice of cutoff radius also plays an important |
| 615 |
< |
role. This connection is further explored below as optimal damping |
| 662 |
< |
coefficients for different choices of $R_\textrm{c}$ are determined |
| 663 |
< |
for {\sc sf} in order to best capture the dielectric behavior. |
| 608 |
> |
``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 |
|
|
| 617 |
< |
\subsubsection{Dielectric Constant}\label{sec:t5peDielectric} |
| 617 |
> |
\subsubsection{Dielectric Constants for TIP5P-E and SSD/RF}\label{sec:t5peDielectric} |
| 618 |
|
|
| 619 |
< |
In the previous, we observed that the choice of damping coefficient |
| 620 |
< |
plays a major role in the calculated dielectric constant. This is not |
| 621 |
< |
too surprising given the results for damping parameter influence on |
| 622 |
< |
the long-time correlated motions of the NaCl crystal in our previous |
| 623 |
< |
study.\cite{Fennell06} The static dielectric constant is calculated |
| 624 |
< |
from the long-time fluctuations of the system's accumulated dipole |
| 625 |
< |
moment (Eq. (\ref{eq:staticDielectric})), so it is going to be quite |
| 626 |
< |
sensitive to the choice of damping parameter. We would like to choose |
| 627 |
< |
optimal damping constants such that any arbitrary choice of cutoff |
| 628 |
< |
radius will properly capture the dielectric behavior of the liquid. |
| 619 |
> |
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 |
|
|
| 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 |
< |
cutoff radius. To calculate the static dielectric constant, we |
| 634 |
< |
performed 5~ns $NPT$ calculations on systems of 512 TIP5P-E water |
| 635 |
< |
molecules and averaged over the converged region (typically the later |
| 636 |
< |
2.5~ns) of equation (\ref{eq:staticDielectric}). The selected cutoff |
| 637 |
< |
radii ranged from 9, 10, 11, and 12~\AA , and the damping parameter |
| 638 |
< |
values ranged from 0.1 to 0.35~\AA$^{-1}$. |
| 633 |
> |
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 |
|
|
| 641 |
|
\begin{table} |
| 642 |
|
\centering |
| 643 |
< |
\caption{TIP5P-E dielectric constant as a function of $\alpha$ and |
| 644 |
< |
$R_\textrm{c}$. The color scale ranges from blue (10) to red (100).} |
| 643 |
> |
\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 |
|
\vspace{6pt} |
| 646 |
< |
\begin{tabular}{ lcccc } |
| 646 |
> |
\begin{tabular}{ lccccccccc } |
| 647 |
|
\toprule |
| 648 |
|
\toprule |
| 649 |
< |
& \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\ |
| 650 |
< |
\cmidrule(lr){2-5} |
| 651 |
< |
$\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 \\ |
| 649 |
> |
& \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 |
|
\midrule |
| 654 |
< |
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 \\ |
| 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} \\ |
| 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 \\ |
| 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 \\ |
| 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 \\ |
| 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 \\ |
| 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 \\ |
| 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} \\ |
| 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 \\ |
| 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} \\ |
| 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 \\ |
| 654 |
> |
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 |
|
\bottomrule |
| 666 |
|
\end{tabular} |
| 667 |
< |
\label{tab:tip5peDielectricMap} |
| 667 |
> |
\label{tab:DielectricMap} |
| 668 |
|
\end{table} |
| 669 |
+ |
|
| 670 |
|
The results of these calculations are displayed in table |
| 671 |
< |
\ref{tab:tip5peDielectricMap}. Coloring of the table is based on the |
| 672 |
< |
numerical values within the cells and was added to better facilitate |
| 673 |
< |
interpretation of the displayed data and highlight trends in the |
| 674 |
< |
dielectric constant behavior. This makes it easy to see that the |
| 675 |
< |
dielectric constant is linear with respect to $\alpha$ and |
| 676 |
< |
$R_\textrm{c}$ in the low to moderate damping regions, and the slope |
| 677 |
< |
is 0.025~\AA$^{-1}\cdot R_\textrm{c}^{-1}$. Another point to note is |
| 678 |
< |
that choosing $\alpha$ and $R_\textrm{c}$ identical to those used with |
| 679 |
< |
the Ewald summation results in the same calculated dielectric |
| 680 |
< |
constant. As an example, in the paper outlining the development of |
| 681 |
< |
TIP5P-E, the real-space cutoff and Ewald coefficient were tethered to |
| 726 |
< |
the system size, and for a 512 molecule system are approximately |
| 727 |
< |
12~\AA\ and 0.25~\AA$^{-1}$ respectively.\cite{Rick04} These |
| 728 |
< |
parameters resulted in a dielectric constant of 92$\pm$14, while with |
| 729 |
< |
{\sc sf} these parameters give a dielectric constant of |
| 671 |
> |
\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 |
|
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 |
< |
constant equal to 63$\pm$1.\cite{Horn04} Calculations using {\sc sf} |
| 686 |
< |
with these parameters and this water model give a dielectric constant |
| 687 |
< |
of 61$\pm$1. Since the dielectric constant is dependent on $\alpha$ |
| 688 |
< |
and $R_\textrm{c}$ with the {\sc sf} technique, it might be |
| 689 |
< |
interesting to investigate the dielectric dependence of the real-space |
| 690 |
< |
Ewald parameters. |
| 685 |
> |
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 |
|
|
| 692 |
< |
Although it is tempting to choose damping parameters equivalent to |
| 693 |
< |
these Ewald examples, the results of our previous work indicate that |
| 694 |
< |
values this high are destructive to both the energetics and |
| 692 |
> |
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 |
|
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 |
| 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 |
< |
behavior for the damped {\sc sf} technique will result in consistent |
| 710 |
> |
behavior for the damped SF technique will result in consistent |
| 711 |
|
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}$ |
| 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 |
< |
(81.5$\pm$1.6).\cite{Mahoney00} As a final note, aside from a slight |
| 719 |
> |
(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 |
|
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 |
< |
To look at the dynamic properties of TIP5P-E when using the {\sc sf} |
| 736 |
> |
To look at the dynamic properties of TIP5P-E when using the SF |
| 737 |
|
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 |
| 799 |
|
\centering |
| 800 |
|
\includegraphics[width=3.5in]{./figures/t5peDynamics.pdf} |
| 801 |
|
\caption{Diffusion constants ({\it upper}) and reorientational time |
| 802 |
< |
constants ({\it lower}) for TIP5P-E using the Ewald sum and {\sc sf} |
| 802 |
> |
constants ({\it lower}) for TIP5P-E using the Ewald sum and SF |
| 803 |
|
technique compared with experiment. Data at temperatures less than |
| 804 |
|
0$^\circ$C were not included in the $\tau_2^y$ plot to allow for |
| 805 |
|
easier comparisons in the more relevant temperature regime.} |
| 808 |
|
Results for the diffusion constants and orientational relaxation times |
| 809 |
|
are shown in figure \ref{fig:t5peDynamics}. From this figure, it is |
| 810 |
|
apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using |
| 811 |
< |
the Ewald sum are reproduced with the {\sc sf} technique. The enhanced |
| 811 |
> |
the Ewald sum are reproduced with the SF technique. The enhanced |
| 812 |
|
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 |
< |
techniques. With the undamped {\sc sf} technique, TIP5P-E tends to |
| 815 |
> |
techniques. With the undamped SF technique, TIP5P-E tends to |
| 816 |
|
diffuse a little faster than with the Ewald sum; however, use of light |
| 817 |
|
to moderate damping results in indistinguishable $D$ values. Though |
| 818 |
< |
not apparent in this figure, {\sc sf} values at the lowest temperature |
| 818 |
> |
not apparent in this figure, SF values at the lowest temperature |
| 819 |
|
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 |
< |
glassy-like phase with the {\sc sf} technique at these lower |
| 822 |
> |
glassy-like phase with the SF technique at these lower |
| 823 |
|
temperatures, though this change seems to be more prominent with the |
| 824 |
< |
{\it undamped} {\sc sf} method, which has stronger local pairwise |
| 824 |
> |
{\it undamped} SF method, which has stronger local pairwise |
| 825 |
|
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 |
< |
experiment fairly well when using the {\sc sf} technique, independent |
| 831 |
> |
experiment fairly well when using the SF technique, independent |
| 832 |
|
of the choice of damping constant. Their are several possible reasons |
| 833 |
|
for this deviation between techniques. The Ewald results were |
| 834 |
< |
calculated using shorter (10ps) trajectories than the {\sc sf} results |
| 835 |
< |
(25ps). A quick calculation from a 10~ps trajectory with {\sc sf} with |
| 834 |
> |
calculated using shorter (10ps) trajectories than the SF results |
| 835 |
> |
(25ps). A quick calculation from a 10~ps trajectory with SF with |
| 836 |
|
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, |
| 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 |
< |
ensemble, while the $NVE$ ensemble was used for {\sc sf} |
| 850 |
> |
ensemble, while the $NVE$ ensemble was used for SF |
| 851 |
|
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 |
| 857 |
|
\subsection{SSD/RF} |
| 858 |
|
|
| 859 |
|
In section \ref{sec:dampingMultipoles}, we described a method for |
| 860 |
< |
applying the damped {\sc sf} technique to systems containing point |
| 860 |
> |
applying the damped SF technique to systems containing point |
| 861 |
|
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 |
| 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 |
< |
consider that the {\sc sf} technique was presented as a pairwise |
| 873 |
> |
consider that the SF technique was presented as a pairwise |
| 874 |
|
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} |
| 921 |
|
SSD/RF can be used effectively with mixed systems, such as dissolved |
| 922 |
|
ions, dissolved organic molecules, or even proteins. |
| 923 |
|
|
| 954 |
– |
\subsubsection{Dielectric Constant} |
| 955 |
– |
|
| 956 |
– |
\begin{table} |
| 957 |
– |
\centering |
| 958 |
– |
\caption{SSD/RF dielectric constant as a function of $\alpha$ and $R_\textrm{c}$. The color scale ranges from blue (10) to red (100).} |
| 959 |
– |
\vspace{6pt} |
| 960 |
– |
\begin{tabular}{ lcccc } |
| 961 |
– |
\toprule |
| 962 |
– |
\toprule |
| 963 |
– |
& \multicolumn{4}{c}{$R_\textrm{c}$ (\AA )} \\ |
| 964 |
– |
\cmidrule(lr){2-5} |
| 965 |
– |
$\alpha$ (\AA$^{-1}$) & 9 & 10 & 11 & 12 \\ |
| 966 |
– |
\midrule |
| 967 |
– |
0.35 & \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 \\ |
| 968 |
– |
& \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} \\ |
| 969 |
– |
0.3 & \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 \\ |
| 970 |
– |
0.275 & \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 \\ |
| 971 |
– |
0.25 & \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 \\ |
| 972 |
– |
0.225 & \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 \\ |
| 973 |
– |
0.2 & \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 \\ |
| 974 |
– |
& \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} \\ |
| 975 |
– |
0.15 & \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 \\ |
| 976 |
– |
& \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} \\ |
| 977 |
– |
0.1 & \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 \\ |
| 978 |
– |
\bottomrule |
| 979 |
– |
\end{tabular} |
| 980 |
– |
\label{tab:ssdrfDielectricMap} |
| 981 |
– |
\end{table} |
| 982 |
– |
In addition to the properties tabulated in table |
| 983 |
– |
\ref{tab:dampedSSDRF}, we mapped out the static dielectric constant of |
| 984 |
– |
SSD/RF as a function of $R_\textrm{c}$ and $\alpha$ in the same |
| 985 |
– |
fashion as in section \ref{sec:t5peDielectric} with TIP5P-E. It is |
| 986 |
– |
apparent from this table that electrostatic damping has a more |
| 987 |
– |
pronounced effect on the dipolar interactions of SSD/RF than the |
| 988 |
– |
monopolar interactions of TIP5P-E. The dielectric constant covers a |
| 989 |
– |
much wider range and has a steeper slope with increasing damping |
| 990 |
– |
parameter. We can use the same trend of $\alpha$ with $R_\textrm{c}$ |
| 991 |
– |
used with TIP5P-E, and for a 12~\AA\ $R_\textrm{c}$, the resulting |
| 992 |
– |
dielectric constant is 82.6$\pm$0.6. This value compares very |
| 993 |
– |
favorably with the experimental value of 78.3.\cite{Malmberg56} This |
| 994 |
– |
is not surprising given that early studies of the SSD model indicate a |
| 995 |
– |
static dielectric constant around 81.\cite{Liu96} |
| 996 |
– |
|
| 924 |
|
\section{Application of Pairwise Electrostatic Corrections: Imaginary Ice} |
| 925 |
|
|
| 926 |
|
In an earlier work, we performed a series of free energy calculations |
| 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 |
< |
observed. Being that the {\sc sf} technique can be used as a simple |
| 932 |
> |
observed. Being that the SF technique can be used as a simple |
| 933 |
|
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 |
< |
discussed ice polymorphs using the {\sc sf} technique with a cutoff |
| 937 |
> |
discussed ice polymorphs using the SF technique with a cutoff |
| 938 |
|
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. |
| 942 |
|
\begin{table} |
| 943 |
|
\centering |
| 944 |
|
\caption{Helmholtz free energies of ice polymorphs at 1~atm and 200~K |
| 945 |
< |
using the damped {\sc sf} electrostatic correction method with a |
| 945 |
> |
using the damped SF electrostatic correction method with a |
| 946 |
|
variety of water models.} |
| 947 |
|
\begin{tabular}{ lccccc } |
| 948 |
|
\toprule |
