| 41 |
|
applies an external pressure to the facets comprising the convex |
| 42 |
|
hull surrounding the objects in the system. Additionally, a Langevin |
| 43 |
|
thermostat is applied to facets of the hull to mimic contact with an |
| 44 |
< |
external heat bath. This new method, the ``Langevin Hull'', |
| 45 |
< |
performs better than traditional affine transform methods for |
| 46 |
< |
systems containing heterogeneous mixtures of materials with |
| 47 |
< |
different compressibilities. It does not suffer from the edge |
| 48 |
< |
effects of boundary potential methods, and allows realistic |
| 49 |
< |
treatment of both external pressure and thermal conductivity to an |
| 50 |
< |
implicit solvents. We apply this method to several different |
| 51 |
< |
systems including bare nano-particles, nano-particles in explicit |
| 52 |
< |
solvent, as well as clusters of liquid water and ice. The predicted |
| 53 |
< |
mechanical and thermal properties of these systems are in good |
| 54 |
< |
agreement with experimental data. |
| 44 |
> |
external heat bath. This new method, the ``Langevin Hull'', performs |
| 45 |
> |
better than traditional affine transform methods for systems |
| 46 |
> |
containing heterogeneous mixtures of materials with different |
| 47 |
> |
compressibilities. It does not suffer from the edge effects of |
| 48 |
> |
boundary potential methods, and allows realistic treatment of both |
| 49 |
> |
external pressure and thermal conductivity to an implicit solvent. |
| 50 |
> |
We apply this method to several different systems including bare |
| 51 |
> |
nanoparticles, nanoparticles in an explicit solvent, as well as |
| 52 |
> |
clusters of liquid water and ice. The predicted mechanical and |
| 53 |
> |
thermal properties of these systems are in good agreement with |
| 54 |
> |
experimental data. |
| 55 |
|
\end{abstract} |
| 56 |
|
|
| 57 |
|
\newpage |
| 75 |
|
system geometry. An affine transform scales both the box lengths as |
| 76 |
|
well as the scaled particle positions (but not the sizes of the |
| 77 |
|
particles). The most common constant pressure methods, including the |
| 78 |
< |
Melchionna modification\cite{melchionna93} to the |
| 79 |
< |
Nos\'e-Hoover-Andersen equations of motion, the Berendsen pressure |
| 80 |
< |
bath, and the Langevin Piston, all utilize coordinate transformation |
| 81 |
< |
to adjust the box volume. |
| 78 |
> |
Melchionna modification\cite{Melchionna1993} to the |
| 79 |
> |
Nos\'e-Hoover-Andersen equations of |
| 80 |
> |
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen |
| 81 |
> |
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin |
| 82 |
> |
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize coordinate |
| 83 |
> |
transformation to adjust the box volume. As long as the material in |
| 84 |
> |
the simulation box is essentially a bulk-like liquid which has a |
| 85 |
> |
relatively uniform compressibility, the standard affine transform |
| 86 |
> |
approach provides an excellent way of adjusting the volume of the |
| 87 |
> |
system and applying pressure directly via the interactions between |
| 88 |
> |
atomic sites. |
| 89 |
|
|
| 90 |
+ |
The problem with this approach becomes apparent when the material |
| 91 |
+ |
being simulated is an inhomogeneous mixture in which portions of the |
| 92 |
+ |
simulation box are incompressible relative to other portions. |
| 93 |
+ |
Examples include simulations of metallic nanoparticles in liquid |
| 94 |
+ |
environments, proteins at interfaces, as well as other multi-phase or |
| 95 |
+ |
interfacial environments. In these cases, the affine transform of |
| 96 |
+ |
atomic coordinates will either cause numerical instability when the |
| 97 |
+ |
sites in the incompressible medium collide with each other, or lead to |
| 98 |
+ |
inefficient sampling of system volumes if the barostat is set slow |
| 99 |
+ |
enough to avoid the instabilities in the incompressible region. |
| 100 |
+ |
|
| 101 |
|
\begin{figure} |
| 102 |
|
\includegraphics[width=\linewidth]{AffineScale2} |
| 103 |
|
\caption{Affine Scaling constant pressure methods use box-length |
| 110 |
|
\label{affineScale} |
| 111 |
|
\end{figure} |
| 112 |
|
|
| 113 |
+ |
One may also wish to avoid affine transform periodic boundary methods |
| 114 |
+ |
to simulate {\it explicitly non-periodic systems} under constant |
| 115 |
+ |
pressure conditions. The use of periodic boxes to enforce a system |
| 116 |
+ |
volume either requires effective solute concentrations that are much |
| 117 |
+ |
higher than desirable, or unreasonable system sizes to avoid this |
| 118 |
+ |
effect. For example, calculations using typical hydration shells |
| 119 |
+ |
solvating a protein under periodic boundary conditions are quite |
| 120 |
+ |
expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL |
| 121 |
+ |
SIMULATIONS] |
| 122 |
|
|
| 123 |
< |
Heterogeneous mixtures of materials with different compressibilities? |
| 123 |
> |
There have been a number of other approaches to explicit |
| 124 |
> |
non-periodicity that focus on constant or nearly-constant {\it volume} |
| 125 |
> |
conditions while maintaining bulk-like behavior. Berkowitz and |
| 126 |
> |
McCammon introduced a stochastic (Langevin) boundary layer inside a |
| 127 |
> |
region of fixed molecules which effectively enforces constant |
| 128 |
> |
temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this |
| 129 |
> |
approach, the stochastic and fixed regions were defined relative to a |
| 130 |
> |
central atom. Brooks and Karplus extended this method to include |
| 131 |
> |
deformable stochastic boundaries.\cite{iii:6312} The stochastic |
| 132 |
> |
boundary approach has been used widely for protein |
| 133 |
> |
simulations. [CITATIONS NEEDED] |
| 134 |
|
|
| 135 |
< |
Explicitly non-periodic systems |
| 135 |
> |
The electrostatic and dispersive behavior near the boundary has long |
| 136 |
> |
been a cause for concern. King and Warshel introduced a surface |
| 137 |
> |
constrained all-atom solvent (SCAAS) which included polarization |
| 138 |
> |
effects of a fixed spherical boundary to mimic bulk-like behavior |
| 139 |
> |
without periodic boundaries.\cite{king:3647} In the SCAAS model, a |
| 140 |
> |
layer of fixed solvent molecules surrounds the solute and any explicit |
| 141 |
> |
solvent, and this in turn is surrounded by a continuum dielectric. |
| 142 |
> |
MORE HERE. WHAT DID THEY FIND? |
| 143 |
|
|
| 144 |
< |
Elastic Bag |
| 144 |
> |
Beglov and Roux developed a boundary model in which the hard sphere |
| 145 |
> |
boundary has a radius that varies with the instantaneous configuration |
| 146 |
> |
of the solute (and solvent) molecules.\cite{beglov:9050} This model |
| 147 |
> |
contains a clear pressure and surface tension contribution to the free |
| 148 |
> |
energy which XXX. |
| 149 |
|
|
| 150 |
< |
Spherical Boundary approaches |
| 150 |
> |
Restraining {\it potentials} introduce repulsive potentials at the |
| 151 |
> |
surface of a sphere or other geometry. The solute and any explicit |
| 152 |
> |
solvent are therefore restrained inside this potential. Often the |
| 153 |
> |
potentials include a weak short-range attraction to maintain the |
| 154 |
> |
correct density at the boundary. Beglov and Roux have also introduced |
| 155 |
> |
a restraining boundary potential which relaxes dynamically depending |
| 156 |
> |
on the solute geometry and the force the explicit system exerts on the |
| 157 |
> |
shell.\cite{Beglov:1995fk} |
| 158 |
|
|
| 159 |
< |
\section{Methodology} |
| 159 |
> |
Recently, Krilov {\it et al.} introduced a flexible boundary model |
| 160 |
> |
that uses a Lennard-Jones potential between the solvent molecules and |
| 161 |
> |
a boundary which is determined dynamically from the position of the |
| 162 |
> |
nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This approach allows |
| 163 |
> |
the confining potential to prevent solvent molecules from migrating |
| 164 |
> |
too far from the solute surface, while providing a weak attractive |
| 165 |
> |
force pulling the solvent molecules towards a fictitious bulk solvent. |
| 166 |
> |
Although this approach is appealing and has physical motivation, |
| 167 |
> |
nanoparticles do not deform far from their original geometries even at |
| 168 |
> |
temperatures which vaporize the nearby solvent. For the systems like |
| 169 |
> |
the one described, the flexible boundary model will be nearly |
| 170 |
> |
identical to a fixed-volume restraining potential. |
| 171 |
|
|
| 172 |
< |
A new method which uses a constant pressure and temperature bath that |
| 173 |
< |
interacts with the objects that are currently at the edge of the |
| 174 |
< |
system. |
| 172 |
> |
The approach of Kohanoff, Caro, and Finnis is the most promising of |
| 173 |
> |
the methods for introducing both constant pressure and temperature |
| 174 |
> |
into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru} |
| 175 |
> |
This method is based on standard Langevin dynamics, but the Brownian |
| 176 |
> |
or random forces are allowed to act only on peripheral atoms and exert |
| 177 |
> |
force in a direction that is inward-facing relative to the facets of a |
| 178 |
> |
closed bounding surface. The statistical distribution of the random |
| 179 |
> |
forces are uniquely tied to the pressure in the external reservoir, so |
| 180 |
> |
the method can be shown to sample the isobaric-isothermal ensemble. |
| 181 |
> |
Kohanoff {\it et al.} used a Delaunay tessellation to generate a |
| 182 |
> |
bounding surface surrounding the outermost atoms in the simulated |
| 183 |
> |
system. This is not the only possible triangulated outer surface, but |
| 184 |
> |
guarantees that all of the random forces point inward towards the |
| 185 |
> |
cluster. |
| 186 |
|
|
| 187 |
< |
Novel features: No a priori geometry is defined, No affine transforms, |
| 188 |
< |
No fictitious particles, No bounding potentials. |
| 187 |
> |
In the following sections, we extend and generalize the approach of |
| 188 |
> |
Kohanoff, Caro, and Finnis. The new method, which we are calling the |
| 189 |
> |
``Langevin Hull'' applies the external pressure, Langevin drag, and |
| 190 |
> |
random forces on the facets of the {\it hull itself} instead of the |
| 191 |
> |
atomic sites comprising the vertices of the hull. This allows us to |
| 192 |
> |
decouple the external pressure contribution from the drag and random |
| 193 |
> |
force. Section \ref{sec:meth} |
| 194 |
|
|
| 195 |
< |
Simulation starts as a collection of atomic locations in 3D (a point |
| 196 |
< |
cloud). |
| 195 |
> |
\section{Methodology} |
| 196 |
> |
\label{sec:meth} |
| 197 |
|
|
| 198 |
< |
Delaunay triangulation finds all facets between coplanar neighbors. |
| 198 |
> |
We have developed a new method which uses a constant pressure and |
| 199 |
> |
temperature bath. This bath interacts only with the objects that are |
| 200 |
> |
currently at the edge of the system. Since the edge is determined |
| 201 |
> |
dynamically as the simulation progresses, no {\it a priori} geometry |
| 202 |
> |
is defined. The pressure and temperature bath interacts {\it |
| 203 |
> |
directly} with the atoms on the edge and not with atoms interior to |
| 204 |
> |
the simulation. This means that there are no affine transforms |
| 205 |
> |
required. There are also no fictitious particles or bounding |
| 206 |
> |
potentials used in this approach. |
| 207 |
|
|
| 208 |
< |
The Convex Hull is the set of facets that have no concave corners at a |
| 209 |
< |
vertex. |
| 208 |
> |
The basics of the method are as follows. The simulation starts as a |
| 209 |
> |
collection of atomic locations in three dimensions (a point cloud). |
| 210 |
> |
Delaunay triangulation is used to find all facets between coplanar |
| 211 |
> |
neighbors. In highly symmetric point clouds, facets can contain many |
| 212 |
> |
atoms, but in all but the most symmetric of cases one might experience |
| 213 |
> |
in a molecular dynamics simulation, the facets are simple triangles in |
| 214 |
> |
3-space that contain exactly three atoms. |
| 215 |
|
|
| 216 |
< |
Molecules on the convex hull are dynamic. As they re-enter the |
| 217 |
< |
cluster, all interactions with the external bath are removed.The |
| 218 |
< |
external bath applies pressure to the facets of the convex hull in |
| 219 |
< |
direct proportion to the area of the facet. Thermal coupling depends on |
| 220 |
< |
the solvent temperature, friction and the size and shape of each |
| 221 |
< |
facet. |
| 127 |
< |
|
| 128 |
< |
\begin{equation} |
| 129 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U |
| 130 |
< |
\end{equation} |
| 216 |
> |
The convex hull is the set of facets that have {\it no concave |
| 217 |
> |
corners} at an atomic site. This eliminates all facets on the |
| 218 |
> |
interior of the point cloud, leaving only those exposed to the |
| 219 |
> |
bath. Sites on the convex hull are dynamic. As molecules re-enter the |
| 220 |
> |
cluster, all interactions between atoms on that molecule and the |
| 221 |
> |
external bath are removed. |
| 222 |
|
|
| 223 |
+ |
For atomic sites in the interior of the point cloud, the equations of |
| 224 |
+ |
motion are simple Newtonian dynamics, |
| 225 |
|
\begin{equation} |
| 226 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext} |
| 226 |
> |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
| 227 |
> |
\label{eq:Newton} |
| 228 |
|
\end{equation} |
| 229 |
+ |
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the |
| 230 |
+ |
instantaneous velocity of site $i$ at time $t$, and $U$ is the total |
| 231 |
+ |
potential energy. For atoms on the exterior of the cluster |
| 232 |
+ |
(i.e. those that occupy one of the vertices of the convex hull), the |
| 233 |
+ |
equation of motion is modified with an external force, ${\mathbf |
| 234 |
+ |
F}_i^{\mathrm ext}$, |
| 235 |
+ |
\begin{equation} |
| 236 |
+ |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
| 237 |
+ |
\end{equation} |
| 238 |
|
|
| 239 |
+ |
The external bath interacts directly with the facets of the convex |
| 240 |
+ |
hull. Since each vertex (or atom) provides one corner of a triangular |
| 241 |
+ |
facet, the force on the facets are divided equally to each vertex. |
| 242 |
+ |
However, each vertex can participate in multiple facets, so the resultant |
| 243 |
+ |
force is a sum over all facets $f$ containing vertex $i$: |
| 244 |
|
\begin{equation} |
| 245 |
|
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\ |
| 246 |
|
} f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf |
| 247 |
|
F}_f^{\mathrm ext} |
| 248 |
|
\end{equation} |
| 249 |
|
|
| 250 |
+ |
The external pressure bath applies a force to the facets of the convex |
| 251 |
+ |
hull in direct proportion to the area of the facet, while the thermal |
| 252 |
+ |
coupling depends on the solvent temperature, friction and the size and |
| 253 |
+ |
shape of each facet. The thermal interactions are expressed as a |
| 254 |
+ |
typical Langevin description of the forces, |
| 255 |
|
\begin{equation} |
| 256 |
|
\begin{array}{rclclcl} |
| 257 |
|
{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\ |
| 258 |
|
& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t) |
| 259 |
|
\end{array} |
| 260 |
|
\end{equation} |
| 261 |
< |
|
| 261 |
> |
Here, $P$ is the external pressure, $A_f$ and $\hat{n}_f$ are the area |
| 262 |
> |
and normal vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is |
| 263 |
> |
the velocity of the facet, |
| 264 |
> |
\begin{equation} |
| 265 |
> |
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
| 266 |
> |
\end{equation} |
| 267 |
> |
and $\Xi_f(t)$ is an approximate ($3 \times 3$) hydrodynamic tensor |
| 268 |
> |
that depends on the geometry and surface area of facet $f$ and the |
| 269 |
> |
viscosity of the fluid (See Appendix A). The hydrodynamic tensor is |
| 270 |
> |
related to the fluctuations of the random force, $\mathbf{R}(t)$, by |
| 271 |
> |
the fluctuation-dissipation theorem, |
| 272 |
|
\begin{eqnarray} |
| 150 |
– |
A_f & = & \text{area of facet}\ f \\ |
| 151 |
– |
\hat{n}_f & = & \text{facet normal} \\ |
| 152 |
– |
P & = & \text{external pressure} |
| 153 |
– |
\end{eqnarray} |
| 154 |
– |
|
| 155 |
– |
\begin{eqnarray} |
| 156 |
– |
{\mathbf v}_f(t) & = & \text{velocity of facet} \\ |
| 157 |
– |
& = & \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i \\ |
| 158 |
– |
\Xi_f(t) & = & \text{is a hydrodynamic tensor that depends} \\ |
| 159 |
– |
& & \text{on the geometry and surface area of} \\ |
| 160 |
– |
& & \text{facet} \ f\ \text{and the viscosity of the fluid.} |
| 161 |
– |
\end{eqnarray} |
| 162 |
– |
|
| 163 |
– |
\begin{eqnarray} |
| 273 |
|
\left< {\mathbf R}_f(t) \right> & = & 0 \\ |
| 274 |
|
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\ |
| 275 |
< |
\Xi_f(t)\delta(t-t^\prime) |
| 275 |
> |
\Xi_f(t)\delta(t-t^\prime). |
| 276 |
> |
\label{eq:randomForce} |
| 277 |
|
\end{eqnarray} |
| 278 |
|
|
| 279 |
< |
Implemented in OpenMD.\cite{Meineke:2005gd,openmd} |
| 279 |
> |
Once the hydrodynamic tensor is known for a given facet (see Appendix |
| 280 |
> |
A) obtaining a stochastic vector that has the properties in |
| 281 |
> |
Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a |
| 282 |
> |
one-time Cholesky decomposition to obtain the square root matrix of |
| 283 |
> |
the resistance tensor, |
| 284 |
> |
\begin{equation} |
| 285 |
> |
\Xi_f = {\bf S} {\bf S}^{T}, |
| 286 |
> |
\label{eq:Cholesky} |
| 287 |
> |
\end{equation} |
| 288 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
| 289 |
> |
vector with the statistics required for the random force can then be |
| 290 |
> |
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which |
| 291 |
> |
has elements chosen from a Gaussian distribution, such that: |
| 292 |
> |
\begin{equation} |
| 293 |
> |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
| 294 |
> |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
| 295 |
> |
\end{equation} |
| 296 |
> |
where $\delta t$ is the timestep in use during the simulation. The |
| 297 |
> |
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
| 298 |
> |
have the correct properties required by Eq. (\ref{eq:randomForce}). |
| 299 |
|
|
| 300 |
+ |
Our treatment of the hydrodynamic tensor must be approximate. $\Xi$ |
| 301 |
+ |
for a triangular plate would normally be treated as a $6 \times 6$ |
| 302 |
+ |
tensor that includes translational and rotational drag as well as |
| 303 |
+ |
translational-rotational coupling. The computation of hydrodynamic |
| 304 |
+ |
tensors for rigid bodies has been detailed |
| 305 |
+ |
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun2008} |
| 306 |
+ |
but the standard approach involving bead approximations would be |
| 307 |
+ |
prohibitively expensive if it were recomputed at each step in a |
| 308 |
+ |
molecular dynamics simulation. |
| 309 |
+ |
|
| 310 |
+ |
We are utilizing an approximate hydrodynamic tensor obtained by first |
| 311 |
+ |
constructing the Oseen tensor for the interaction of the centroid of |
| 312 |
+ |
the facet ($f$) with each of the subfacets $j$, |
| 313 |
+ |
\begin{equation} |
| 314 |
+ |
T_{jf}=\frac{A_j}{8\pi\eta R_{jf}}\left(I + |
| 315 |
+ |
\frac{\mathbf{R}_{jf}\mathbf{R}_{jf}^T}{R_{jf}^2}\right) |
| 316 |
+ |
\end{equation} |
| 317 |
+ |
Here, $A_j$ is the area of subfacet $j$ which is a triangle containing |
| 318 |
+ |
two of the vertices of the facet along with the centroid. |
| 319 |
+ |
$\mathbf{R}_{jf}$ is the vector between the centroid of facet $f$ and |
| 320 |
+ |
the centroid of sub-facet $j$, and $I$ is the ($3 \times 3$) identity |
| 321 |
+ |
matrix. $\eta$ is the viscosity of the external bath. |
| 322 |
+ |
|
| 323 |
+ |
\begin{figure} |
| 324 |
+ |
\includegraphics[width=\linewidth]{hydro} |
| 325 |
+ |
\caption{The hydrodynamic tensor $\Xi$ for a facet comprising sites $i$, |
| 326 |
+ |
$j$, and $k$ is constructed using Oseen tensor contributions |
| 327 |
+ |
between the centoid of the facet $f$ and each of the sub-facets |
| 328 |
+ |
($i,f,j$), ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets |
| 329 |
+ |
are located at $1$, $2$, and $3$, and the area of each sub-facet is |
| 330 |
+ |
easily computed using half the cross product of two of the edges.} |
| 331 |
+ |
\label{hydro} |
| 332 |
+ |
\end{figure} |
| 333 |
+ |
|
| 334 |
+ |
The Oseen tensors for each of the sub-facets are summed, and the |
| 335 |
+ |
resulting matrix is inverted to give a $3 \times 3$ hydrodynamic |
| 336 |
+ |
tensor for translations of the triangular plate, |
| 337 |
+ |
\begin{equation} |
| 338 |
+ |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
| 339 |
+ |
\end{equation} |
| 340 |
+ |
We have implemented this method by extending the Langevin dynamics |
| 341 |
+ |
integrator in our group code, OpenMD.\cite{Meineke2005,openmd} There |
| 342 |
+ |
is a moderate penalty for computing the convex hull at each step in |
| 343 |
+ |
the molecular dynamics simulation (HOW MUCH?), but the convex hull is |
| 344 |
+ |
remarkably easy to parallelize on distributed memory machines (see |
| 345 |
+ |
Appendix B). |
| 346 |
+ |
|
| 347 |
|
\section{Tests \& Applications} |
| 348 |
+ |
\label{sec:tests} |
| 349 |
|
|
| 350 |
|
\subsection{Bulk modulus of gold nanoparticles} |
| 351 |
|
|
| 384 |
|
\label{compWater} |
| 385 |
|
\end{figure} |
| 386 |
|
|
| 387 |
< |
We initially used the classic compressibility formula |
| 387 |
> |
The volume of a three-dimensional point cloud is not an obvious property to calculate. In order to calculate the isothermal compressibility we adapted the classic compressibility formula so that the compressibility could be calculated using information about the local density instead of the total volume of the convex hull. |
| 388 |
|
|
| 389 |
|
\begin{equation} |
| 390 |
|
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T} |
| 391 |
|
\end{equation} |
| 392 |
|
|
| 216 |
– |
to calculate the the isothermal compressibility at each target pressure. These calculations yielded compressibility values that were dramatically higher than both previous simulations and experiment. The particular compressibility expression used requires the calculation of both a volume and pressure differential, thereby stipulating that the data from at least two simulations at different pressures must be used to calculate the isothermal compressibility at one pressure. |
| 393 |
|
|
| 394 |
< |
Per the fluctuation dissipation theorem \cite{Debendedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure |
| 394 |
> |
Assuming a uniform density, we can use the relationship $\rho = \frac{N}{V}$ to rewrite the isothermal compressibility formula as |
| 395 |
|
|
| 396 |
|
\begin{equation} |
| 397 |
< |
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T} |
| 397 |
> |
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T} |
| 398 |
|
\end{equation} |
| 399 |
|
|
| 400 |
< |
Thus, the compressibility of each simulation run can be calculated entirely independently from all other trajectories. However, the resulting compressibilities were still as much as an order of magnitude larger than the reference values. The effect was particularly pronounced at the low end of the pressure range. At ambient temperature and low pressures, there exists an equilibrium between vapor and liquid phases. Vapor molecules are naturally more diffuse around the exterior of the cluster, causing artificially large cluster volumes. Any compressibility calculation that relies on the hull volume will suffer these effects. |
| 400 |
> |
Isothermal compressibility values calculated using this modified expression are in good agreement with the reference values throughout the 1 - 1000 atm pressure regime. Regardless of the difficulty in obtaining accurate hull volumes at low temperature and pressures, the Langevin Hull NPT method provides reasonable isothermal compressibility values for water through a large range of pressures. |
| 401 |
|
|
| 402 |
< |
In order to calculate the isothermal compressibility without being hindered by hull volume issues, we adapted the classic compressibility formula so that the compressibility could be calculated using information about the local density instead of the volume of the convex hull. We calculated the $g_{OO}(r)$ for a 1 nanosecond simulation of a cluster of 1372 SPC/E water molecules and spherically integrated the function over the bounds 0 to $r'$. In all cases, the value of $r'$ was 17.26216 $\AA$. The value of the total integral between these bounds is essentially the number (N) of molecules within volume $\frac{4}{3}\pi r'^{3}$ at a given pressure. To yield an actual molecule count, N must be scaled by an ideal density. However, even in the absence of an ideal density, we can use the relationship $\rho = \frac{N}{V}$ to rewrite the isothermal compressibility formula as |
| 402 |
> |
We initially used the classic compressibility formula to calculate the the isothermal compressibility at each target pressure. These calculations yielded compressibility values that were dramatically higher than both previous simulations and experiment. The particular compressibility expression used requires the calculation of both a volume and pressure differential, thereby stipulating that the data from at least two simulations at different pressures must be used to calculate the isothermal compressibility at one pressure. |
| 403 |
|
|
| 404 |
+ |
Per the fluctuation dissipation theorem \cite{Debenedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure |
| 405 |
+ |
|
| 406 |
|
\begin{equation} |
| 407 |
< |
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T} |
| 407 |
> |
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T} |
| 408 |
|
\end{equation} |
| 409 |
|
|
| 410 |
< |
Isothermal compressibility values calculated using this modified expression are in good agreement with the reference values throughout the 1 - 1000 atm pressure regime. Regardless of the difficulty in obtaining accurate hull volumes at low temperature and pressures, the Langevin Hull NPT method provides reasonable isothermal compressibility values for water through a large range of pressures. |
| 410 |
> |
Thus, the compressibility of each simulation run can be calculated entirely independently from all other trajectories. However, the resulting compressibilities were still as much as an order of magnitude larger than the reference values. The effect was particularly pronounced at the low end of the pressure range. At ambient temperature and low pressures, there exists an equilibrium between vapor and liquid phases. Vapor molecules are naturally more diffuse around the exterior of the cluster, causing artificially large cluster volumes. Any compressibility calculation that relies on the hull volume will suffer these effects. |
| 411 |
|
|
| 412 |
+ |
|
| 413 |
|
\subsection{Molecular orientation distribution at cluster boundary} |
| 414 |
|
|
| 415 |
|
In order for non-periodic boundary conditions to be widely applicable, they must be constructed in such a way that they allow a finite, usually small, simulated system to replicate the properties of an infinite bulk system. Naturally, this requirement has spawned many methods for inserting boundaries into simulated systems [REF... ?]. Of particular interest to our characterization of the Langevin Hull is the orientation of water molecules included in the geometric hull. Ideally, all molecules in the cluster will have the same orientational distribution as bulk water. |
| 454 |
|
|
| 455 |
|
\section{Appendix A: Hydrodynamic tensor for triangular facets} |
| 456 |
|
|
| 278 |
– |
\begin{figure} |
| 279 |
– |
\includegraphics[width=\linewidth]{hydro} |
| 280 |
– |
\caption{Hydro} |
| 281 |
– |
\label{hydro} |
| 282 |
– |
\end{figure} |
| 283 |
– |
|
| 284 |
– |
\begin{equation} |
| 285 |
– |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1} |
| 286 |
– |
\end{equation} |
| 287 |
– |
|
| 288 |
– |
\begin{equation} |
| 289 |
– |
T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I + |
| 290 |
– |
\frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right) |
| 291 |
– |
\end{equation} |
| 292 |
– |
|
| 457 |
|
\section{Appendix B: Computing Convex Hulls on Parallel Computers} |
| 458 |
|
|
| 459 |
|
\section{Acknowledgments} |
