| 17 |
|
\setlength{\abovecaptionskip}{20 pt} |
| 18 |
|
\setlength{\belowcaptionskip}{30 pt} |
| 19 |
|
|
| 20 |
< |
\bibpunct{[}{]}{,}{s}{}{;} |
| 20 |
> |
\bibpunct{}{}{,}{s}{}{;} |
| 21 |
|
\bibliographystyle{achemso} |
| 22 |
|
|
| 23 |
|
\begin{document} |
| 42 |
|
hull surrounding the system. A Langevin thermostat is also applied |
| 43 |
|
to the facets to mimic contact with an external heat bath. This new |
| 44 |
|
method, the ``Langevin Hull'', can handle heterogeneous mixtures of |
| 45 |
< |
materials with different compressibilities. These are systems that |
| 46 |
< |
are problematic for traditional affine transform methods. The |
| 47 |
< |
Langevin Hull does not suffer from the edge effects of boundary |
| 48 |
< |
potential methods, and allows realistic treatment of both external |
| 49 |
< |
pressure and thermal conductivity due to the presence of an implicit |
| 50 |
< |
solvent. We apply this method to several different systems |
| 51 |
< |
including bare metal nanoparticles, nanoparticles in an explicit |
| 52 |
< |
solvent, as well as clusters of liquid water. The predicted |
| 53 |
< |
mechanical properties of these systems are in good agreement with |
| 54 |
< |
experimental data and previous simulation work. |
| 45 |
> |
materials with different compressibilities. These systems are |
| 46 |
> |
problematic for traditional affine transform methods. The Langevin |
| 47 |
> |
Hull does not suffer from the edge effects of boundary potential |
| 48 |
> |
methods, and allows realistic treatment of both external pressure |
| 49 |
> |
and thermal conductivity due to the presence of an implicit solvent. |
| 50 |
> |
We apply this method to several different systems including bare |
| 51 |
> |
metal nanoparticles, nanoparticles in an explicit solvent, as well |
| 52 |
> |
as clusters of liquid water. The predicted mechanical properties of |
| 53 |
> |
these systems are in good agreement with experimental data and |
| 54 |
> |
previous simulation work. |
| 55 |
|
\end{abstract} |
| 56 |
|
|
| 57 |
|
\newpage |
| 117 |
|
higher than desirable, or unreasonable system sizes to avoid this |
| 118 |
|
effect. For example, calculations using typical hydration boxes |
| 119 |
|
solvating a protein under periodic boundary conditions are quite |
| 120 |
< |
expensive. A 62 $\AA^3$ box of water solvating a moderately small |
| 120 |
> |
expensive. A 62 \AA$^3$ box of water solvating a moderately small |
| 121 |
|
protein like hen egg white lysozyme (PDB code: 1LYZ) yields an |
| 122 |
|
effective protein concentration of 100 mg/mL.\cite{Asthagiri20053300} |
| 123 |
|
|
| 124 |
< |
Typically protein concentrations in the cell are on the order of |
| 125 |
< |
160-310 mg/ml,\cite{Brown1991195} and the factor of 20 difference |
| 126 |
< |
between simulations and the cellular environment may have significant |
| 127 |
< |
effects on the structure and dynamics of simulated protein structures. |
| 124 |
> |
{\it Total} protein concentrations in the cell are typically on the |
| 125 |
> |
order of 160-310 mg/ml,\cite{Brown1991195} and individual proteins |
| 126 |
> |
have concentrations orders of magnitude lower than this in the |
| 127 |
> |
cellular environment. The effective concentrations of single proteins |
| 128 |
> |
in simulations may have significant effects on the structure and |
| 129 |
> |
dynamics of simulated systems. |
| 130 |
|
|
| 129 |
– |
|
| 131 |
|
\subsection*{Boundary Methods} |
| 132 |
|
There have been a number of approaches to handle simulations of |
| 133 |
|
explicitly non-periodic systems that focus on constant or |
| 139 |
|
fixed regions were defined relative to a central atom. Brooks and |
| 140 |
|
Karplus extended this method to include deformable stochastic |
| 141 |
|
boundaries.\cite{iii:6312} The stochastic boundary approach has been |
| 142 |
< |
used widely for protein simulations. [CITATIONS NEEDED] |
| 142 |
> |
used widely for protein simulations. |
| 143 |
|
|
| 144 |
|
The electrostatic and dispersive behavior near the boundary has long |
| 145 |
|
been a cause for concern when performing simulations of explicitly |
| 150 |
|
simulated clusters of TIPS2 water surrounded by a hydrophobic bounding |
| 151 |
|
potential. The spherical hydrophobic boundary induced dangling |
| 152 |
|
hydrogen bonds at the surface that propagated deep into the cluster, |
| 153 |
< |
affecting most of molecules in the simulation. This result echoes an |
| 154 |
< |
earlier study which showed that an extended planar hydrophobic surface |
| 155 |
< |
caused orientational preference at the surface which extended |
| 156 |
< |
relatively deep (7 \r{A}) into the liquid simulation |
| 157 |
< |
cell.\cite{Lee1984} The surface constrained all-atom solvent (SCAAS) |
| 158 |
< |
model \cite{King1989} improved upon its SCSSD predecessor. The SCAAS |
| 159 |
< |
model utilizes a polarization constraint which is applied to the |
| 160 |
< |
surface molecules to maintain bulk-like structure at the cluster |
| 161 |
< |
surface. A radial constraint is used to maintain the desired bulk |
| 162 |
< |
density of the liquid. Both constraint forces are applied only to a |
| 163 |
< |
pre-determined number of the outermost molecules. |
| 153 |
> |
affecting most of the molecules in the simulation. This result echoes |
| 154 |
> |
an earlier study which showed that an extended planar hydrophobic |
| 155 |
> |
surface caused orientational preferences at the surface which extended |
| 156 |
> |
relatively deep (7 \AA) into the liquid simulation cell.\cite{Lee1984} |
| 157 |
> |
The surface constrained all-atom solvent (SCAAS) model \cite{King1989} |
| 158 |
> |
improved upon its SCSSD predecessor. The SCAAS model utilizes a |
| 159 |
> |
polarization constraint which is applied to the surface molecules to |
| 160 |
> |
maintain bulk-like structure at the cluster surface. A radial |
| 161 |
> |
constraint is used to maintain the desired bulk density of the |
| 162 |
> |
liquid. Both constraint forces are applied only to a pre-determined |
| 163 |
> |
number of the outermost molecules. |
| 164 |
|
|
| 165 |
|
Beglov and Roux have developed a boundary model in which the hard |
| 166 |
|
sphere boundary has a radius that varies with the instantaneous |
| 167 |
|
configuration of the solute (and solvent) molecules.\cite{beglov:9050} |
| 168 |
|
This model contains a clear pressure and surface tension contribution |
| 169 |
< |
to the free energy which XXX. |
| 169 |
> |
to the free energy. |
| 170 |
|
|
| 171 |
|
\subsection*{Restraining Potentials} |
| 172 |
|
Restraining {\it potentials} introduce repulsive potentials at the |
| 181 |
|
Recently, Krilov {\it et al.} introduced a {\it flexible} boundary |
| 182 |
|
model that uses a Lennard-Jones potential between the solvent |
| 183 |
|
molecules and a boundary which is determined dynamically from the |
| 184 |
< |
position of the nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This |
| 184 |
> |
position of the nearest solute atom.\cite{LiY._jp046852t,Zhu:2008fk} This |
| 185 |
|
approach allows the confining potential to prevent solvent molecules |
| 186 |
|
from migrating too far from the solute surface, while providing a weak |
| 187 |
|
attractive force pulling the solvent molecules towards a fictitious |
| 197 |
|
into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru} |
| 198 |
|
This method is based on standard Langevin dynamics, but the Brownian |
| 199 |
|
or random forces are allowed to act only on peripheral atoms and exert |
| 200 |
< |
force in a direction that is inward-facing relative to the facets of a |
| 201 |
< |
closed bounding surface. The statistical distribution of the random |
| 200 |
> |
forces in a direction that is inward-facing relative to the facets of |
| 201 |
> |
a closed bounding surface. The statistical distribution of the random |
| 202 |
|
forces are uniquely tied to the pressure in the external reservoir, so |
| 203 |
|
the method can be shown to sample the isobaric-isothermal ensemble. |
| 204 |
|
Kohanoff {\it et al.} used a Delaunay tessellation to generate a |
| 222 |
|
\label{sec:meth} |
| 223 |
|
|
| 224 |
|
The Langevin Hull uses an external bath at a fixed constant pressure |
| 225 |
< |
($P$) and temperature ($T$). This bath interacts only with the |
| 226 |
< |
objects on the exterior hull of the system. Defining the hull of the |
| 227 |
< |
simulation is done in a manner similar to the approach of Kohanoff, |
| 228 |
< |
Caro and Finnis.\cite{Kohanoff:2005qm} That is, any instantaneous |
| 229 |
< |
configuration of the atoms in the system is considered as a point |
| 230 |
< |
cloud in three dimensional space. Delaunay triangulation is used to |
| 231 |
< |
find all facets between coplanar |
| 232 |
< |
neighbors.\cite{delaunay,springerlink:10.1007/BF00977785} In highly |
| 225 |
> |
($P$) and temperature ($T$) with an effective solvent viscosity |
| 226 |
> |
($\eta$). This bath interacts only with the objects on the exterior |
| 227 |
> |
hull of the system. Defining the hull of the atoms in a simulation is |
| 228 |
> |
done in a manner similar to the approach of Kohanoff, Caro and |
| 229 |
> |
Finnis.\cite{Kohanoff:2005qm} That is, any instantaneous configuration |
| 230 |
> |
of the atoms in the system is considered as a point cloud in three |
| 231 |
> |
dimensional space. Delaunay triangulation is used to find all facets |
| 232 |
> |
between coplanar |
| 233 |
> |
neighbors.\cite{delaunay,springerlink:10.1007/BF00977785} In highly |
| 234 |
|
symmetric point clouds, facets can contain many atoms, but in all but |
| 235 |
< |
the most symmetric of cases the facets are simple triangles in 3-space |
| 236 |
< |
that contain exactly three atoms. |
| 235 |
> |
the most symmetric of cases, the facets are simple triangles in |
| 236 |
> |
3-space which contain exactly three atoms. |
| 237 |
|
|
| 238 |
|
The convex hull is the set of facets that have {\it no concave |
| 239 |
|
corners} at an atomic site.\cite{Barber96,EDELSBRUNNER:1994oq} This |
| 247 |
|
simulation. |
| 248 |
|
|
| 249 |
|
\begin{figure} |
| 250 |
< |
\includegraphics[width=\linewidth]{hullSample} |
| 250 |
> |
\includegraphics[width=\linewidth]{solvatedNano} |
| 251 |
|
\caption{The external temperature and pressure bath interacts only |
| 252 |
|
with those atoms on the convex hull (grey surface). The hull is |
| 253 |
|
computed dynamically at each time step, and molecules can move |
| 254 |
< |
between the interior (Newtonian) region and the Langevin hull.} |
| 254 |
> |
between the interior (Newtonian) region and the Langevin Hull.} |
| 255 |
|
\label{fig:hullSample} |
| 256 |
|
\end{figure} |
| 257 |
|
|
| 266 |
|
potential energy. For atoms on the exterior of the cluster |
| 267 |
|
(i.e. those that occupy one of the vertices of the convex hull), the |
| 268 |
|
equation of motion is modified with an external force, ${\mathbf |
| 269 |
< |
F}_i^{\mathrm ext}$, |
| 269 |
> |
F}_i^{\mathrm ext}$: |
| 270 |
|
\begin{equation} |
| 271 |
|
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
| 272 |
|
\end{equation} |
| 302 |
|
\end{equation} |
| 303 |
|
and $\Xi_f(t)$ is an approximate ($3 \times 3$) resistance tensor that |
| 304 |
|
depends on the geometry and surface area of facet $f$ and the |
| 305 |
< |
viscosity of the fluid. The resistance tensor is related to the |
| 305 |
> |
viscosity of the bath. The resistance tensor is related to the |
| 306 |
|
fluctuations of the random force, $\mathbf{R}(t)$, by the |
| 307 |
|
fluctuation-dissipation theorem, |
| 308 |
|
\begin{eqnarray} |
| 332 |
|
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
| 333 |
|
have the correct properties required by Eq. (\ref{eq:randomForce}). |
| 334 |
|
|
| 335 |
< |
Our treatment of the resistance tensor is approximate. $\Xi$ for a |
| 335 |
> |
Our treatment of the resistance tensor is approximate. $\Xi_f$ for a |
| 336 |
|
rigid triangular plate would normally be treated as a $6 \times 6$ |
| 337 |
|
tensor that includes translational and rotational drag as well as |
| 338 |
|
translational-rotational coupling. The computation of resistance |
| 378 |
|
configurations, so this appears to be a reasonably good approximation. |
| 379 |
|
|
| 380 |
|
We have implemented this method by extending the Langevin dynamics |
| 381 |
< |
integrator in our code, OpenMD.\cite{Meineke2005,openmd} At each |
| 381 |
> |
integrator in our code, OpenMD.\cite{Meineke2005,open_md} At each |
| 382 |
|
molecular dynamics time step, the following process is carried out: |
| 383 |
|
\begin{enumerate} |
| 384 |
|
\item The standard inter-atomic forces ($\nabla_iU$) are computed. |
| 387 |
|
\item The convex hull is computed and facets are identified. |
| 388 |
|
\item For each facet: |
| 389 |
|
\begin{itemize} |
| 390 |
< |
\item[a.] The force from the pressure bath ($-PA_f\hat{n}_f$) is |
| 390 |
> |
\item[a.] The force from the pressure bath ($-\hat{n}_fPA_f$) is |
| 391 |
|
computed. |
| 392 |
|
\item[b.] The resistance tensor ($\Xi_f(t)$) is computed using the |
| 393 |
|
viscosity ($\eta$) of the bath. |
| 400 |
|
\item Atomic positions and velocities are propagated. |
| 401 |
|
\end{enumerate} |
| 402 |
|
The Delaunay triangulation and computation of the convex hull are done |
| 403 |
< |
using calls to the qhull library.\cite{Qhull} There is a minimal |
| 403 |
> |
using calls to the qhull library.\cite{Q_hull} There is a minimal |
| 404 |
|
penalty for computing the convex hull and resistance tensors at each |
| 405 |
|
step in the molecular dynamics simulation (roughly 0.02 $\times$ cost |
| 406 |
|
of a single force evaluation), and the convex hull is remarkably easy |
| 412 |
|
To test the new method, we have carried out simulations using the |
| 413 |
|
Langevin Hull on: 1) a crystalline system (gold nanoparticles), 2) a |
| 414 |
|
liquid droplet (SPC/E water),\cite{Berendsen1987} and 3) a |
| 415 |
< |
heterogeneous mixture (gold nanoparticles in a water droplet). In each |
| 414 |
< |
case, we have computed properties that depend on the external applied |
| 415 |
< |
pressure. Of particular interest for the single-phase systems is the |
| 416 |
< |
isothermal compressibility, |
| 415 |
> |
heterogeneous mixture (gold nanoparticles in an SPC/E water droplet). In each case, we have computed properties that depend on the external applied pressure. Of particular interest for the single-phase systems is the isothermal compressibility, |
| 416 |
|
\begin{equation} |
| 417 |
|
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right |
| 418 |
|
)_{T}. |
| 421 |
|
|
| 422 |
|
One problem with eliminating periodic boundary conditions and |
| 423 |
|
simulation boxes is that the volume of a three-dimensional point cloud |
| 424 |
< |
is not well-defined. In order to compute the compressibility of a |
| 424 |
> |
is not well-defined. In order to compute the compressibility of a |
| 425 |
|
bulk material, we make an assumption that the number density, $\rho = |
| 426 |
< |
\frac{N}{V}$, is uniform within some region of the point cloud. The |
| 426 |
> |
\frac{N}{V}$, is uniform within some region of the point cloud. The |
| 427 |
|
compressibility can then be expressed in terms of the average number |
| 428 |
|
of particles in that region, |
| 429 |
|
\begin{equation} |
| 430 |
|
\kappa_{T} = -\frac{1}{N} \left ( \frac{\partial N}{\partial P} \right |
| 431 |
< |
)_{T} |
| 431 |
> |
)_{T}. |
| 432 |
|
\label{eq:BMN} |
| 433 |
|
\end{equation} |
| 434 |
< |
The region we used is a spherical volume of 10 \AA\ radius centered in |
| 435 |
< |
the middle of the cluster. $N$ is the average number of molecules |
| 434 |
> |
The region we used is a spherical volume of 20 \AA\ radius centered in |
| 435 |
> |
the middle of the cluster with a roughly 25 \AA\ radius. $N$ is the average number of molecules |
| 436 |
|
found within this region throughout a given simulation. The geometry |
| 437 |
< |
and size of the region is arbitrary, and any bulk-like portion of the |
| 438 |
< |
cluster can be used to compute the compressibility. |
| 437 |
> |
of the region is arbitrary, and any bulk-like portion of the |
| 438 |
> |
cluster can be used to compute the compressibility. |
| 439 |
|
|
| 440 |
|
One might assume that the volume of the convex hull could simply be |
| 441 |
|
taken as the system volume $V$ in the compressibility expression |
| 453 |
|
Spohr potential was adopted in depicting the interaction between metal |
| 454 |
|
atoms and the SPC/E water molecules.\cite{ISI:000167766600035} |
| 455 |
|
|
| 456 |
< |
\subsection{Compressibility of gold nanoparticles} |
| 456 |
> |
\subsection{Bulk Modulus of gold nanoparticles} |
| 457 |
|
|
| 458 |
|
The compressibility (and its inverse, the bulk modulus) is well-known |
| 459 |
|
for gold, and is captured well by the embedded atom method |
| 460 |
< |
(EAM)~\cite{PhysRevB.33.7983} potential |
| 461 |
< |
and related multi-body force fields. In particular, the quantum |
| 462 |
< |
Sutton-Chen potential gets nearly quantitative agreement with the |
| 463 |
< |
experimental bulk modulus values, and makes a good first test of how |
| 464 |
< |
the Langevin Hull will perform at large applied pressures. |
| 460 |
> |
(EAM)~\cite{PhysRevB.33.7983} potential and related multi-body force |
| 461 |
> |
fields. In particular, the quantum Sutton-Chen potential gets nearly |
| 462 |
> |
quantitative agreement with the experimental bulk modulus values, and |
| 463 |
> |
makes a good first test of how the Langevin Hull will perform at large |
| 464 |
> |
applied pressures. |
| 465 |
|
|
| 466 |
|
The Sutton-Chen (SC) potentials are based on a model of a metal which |
| 467 |
|
treats the nuclei and core electrons as pseudo-atoms embedded in the |
| 468 |
|
electron density due to the valence electrons on all of the other |
| 469 |
< |
atoms in the system.\cite{Chen90} The SC potential has a simple form that closely |
| 470 |
< |
resembles the Lennard Jones potential, |
| 469 |
> |
atoms in the system.\cite{Chen90} The SC potential has a simple form |
| 470 |
> |
that closely resembles the Lennard Jones potential, |
| 471 |
|
\begin{equation} |
| 472 |
|
\label{eq:SCP1} |
| 473 |
|
U_{tot}=\sum _{i}\left[ \frac{1}{2}\sum _{j\neq i}D_{ij}V^{pair}_{ij}(r_{ij})-c_{i}D_{ii}\sqrt{\rho_{i}}\right] , |
| 481 |
|
interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in |
| 482 |
|
Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models |
| 483 |
|
the interactions between the valence electrons and the cores of the |
| 484 |
< |
pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy |
| 484 |
> |
pseudo-atoms. $D_{ij}$ and $D_{ii}$ set the appropriate overall energy |
| 485 |
|
scale, $c_i$ scales the attractive portion of the potential relative |
| 486 |
|
to the repulsive interaction and $\alpha_{ij}$ is a length parameter |
| 487 |
|
that assures a dimensionless form for $\rho$. These parameters are |
| 489 |
|
energy, and elastic moduli for FCC transition metals. The quantum |
| 490 |
|
Sutton-Chen (QSC) formulation matches these properties while including |
| 491 |
|
zero-point quantum corrections for different transition |
| 492 |
< |
metals.\cite{PhysRevB.59.3527} |
| 492 |
> |
metals.\cite{PhysRevB.59.3527,QSC2} |
| 493 |
|
|
| 494 |
|
In bulk gold, the experimentally-measured value for the bulk modulus |
| 495 |
|
is 180.32 GPa, while previous calculations on the QSC potential in |
| 496 |
< |
periodic-boundary simulations of the bulk have yielded values of |
| 497 |
< |
175.53 GPa.\cite{XXX} Using the same force field, we have performed a |
| 498 |
< |
series of relatively short (200 ps) simulations on 40 \r{A} radius |
| 499 |
< |
nanoparticles under the Langevin Hull at a variety of applied |
| 500 |
< |
pressures ranging from 0 GPa to XXX. We obtain a value of 177.547 GPa |
| 501 |
< |
for the bulk modulus for gold using this echnique. |
| 496 |
> |
periodic-boundary simulations of the bulk crystal have yielded values |
| 497 |
> |
of 175.53 GPa.\cite{QSC2} Using the same force field, we have |
| 498 |
> |
performed a series of 1 ns simulations on gold nanoparticles of three |
| 499 |
> |
different radii: 20 \AA~ (1985 atoms), 30 \AA~ (6699 atoms), and 40 |
| 500 |
> |
\AA~ (15707 atoms) utilizing the Langevin Hull at a variety of applied |
| 501 |
> |
pressures ranging from 0 -- 10 GPa. For the 40 \AA~ radius |
| 502 |
> |
nanoparticle we obtain a value of 177.55 GPa for the bulk modulus of |
| 503 |
> |
gold, in close agreement with both previous simulations and the |
| 504 |
> |
experimental bulk modulus reported for gold single |
| 505 |
> |
crystals.\cite{Collard1991} The smaller gold nanoparticles (30 and 20 |
| 506 |
> |
\AA~ radii) have calculated bulk moduli of 215.58 and 208.86 GPa, |
| 507 |
> |
respectively, indicating that smaller nanoparticles are somewhat |
| 508 |
> |
stiffer (less compressible) than the larger nanoparticles. This |
| 509 |
> |
stiffening of the small nanoparticles may be related to their high |
| 510 |
> |
degree of surface curvature, resulting in a lower coordination number |
| 511 |
> |
of surface atoms relative to the the surface atoms in the 40 \AA~ |
| 512 |
> |
radius particle. |
| 513 |
|
|
| 514 |
+ |
We obtain a gold lattice constant of 4.051 \AA~ using the Langevin |
| 515 |
+ |
Hull at 1 atm, close to the experimentally-determined value for bulk |
| 516 |
+ |
gold and the value for gold simulated using the QSC potential and |
| 517 |
+ |
periodic boundary conditions (4.079 \AA~ and 4.088\AA~, |
| 518 |
+ |
respectively).\cite{QSC2} The slightly smaller calculated lattice |
| 519 |
+ |
constant is most likely due to the presence of surface tension in the |
| 520 |
+ |
non-periodic Langevin Hull cluster, an effect absent from a bulk |
| 521 |
+ |
simulation. The specific heat of a 40 \AA~ gold nanoparticle under the |
| 522 |
+ |
Langevin Hull at 1 atm is 24.914 $\mathrm {\frac{J}{mol \, K}}$, which |
| 523 |
+ |
compares very well with the experimental value of 25.42 $\mathrm |
| 524 |
+ |
{\frac{J}{mol \, K}}$. |
| 525 |
+ |
|
| 526 |
|
\begin{figure} |
| 527 |
|
\includegraphics[width=\linewidth]{stacked} |
| 528 |
|
\caption{The response of the internal pressure and temperature of gold |
| 530 |
|
($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting |
| 531 |
|
from initial conditions that were far from the bath pressure and |
| 532 |
|
temperature. The pressure response is rapid (after the breathing mode oscillations in the nanoparticle die out), and the rate of thermal equilibration depends on both exposed surface area (top panel) and the viscosity of the bath (middle panel).} |
| 533 |
< |
\label{pressureResponse} |
| 533 |
> |
\label{fig:pressureResponse} |
| 534 |
|
\end{figure} |
| 535 |
|
|
| 536 |
< |
\begin{equation} |
| 537 |
< |
\kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial |
| 538 |
< |
P}\right) |
| 539 |
< |
\end{equation} |
| 536 |
> |
We note that the Langevin Hull produces rapidly-converging behavior |
| 537 |
> |
for structures that are started far from equilibrium. In |
| 538 |
> |
Fig. \ref{fig:pressureResponse} we show how the pressure and |
| 539 |
> |
temperature respond to the Langevin Hull for nanoparticles that were |
| 540 |
> |
initialized far from the target pressure and temperature. As |
| 541 |
> |
expected, the rate at which thermal equilibrium is achieved depends on |
| 542 |
> |
the total surface area of the cluster exposed to the bath as well as |
| 543 |
> |
the bath viscosity. Pressure that is applied suddenly to a cluster |
| 544 |
> |
can excite breathing vibrations, but these rapidly damp out (on time |
| 545 |
> |
scales of 30 -- 50 ps). |
| 546 |
|
|
| 547 |
|
\subsection{Compressibility of SPC/E water clusters} |
| 548 |
|
|
| 551 |
|
ensembles) have yielded values for the isothermal compressibility that |
| 552 |
|
agree well with experiment.\cite{Fine1973} The results of two |
| 553 |
|
different approaches for computing the isothermal compressibility from |
| 554 |
< |
Langevin Hull simulations for pressures between 1 and 6500 atm are |
| 554 |
> |
Langevin Hull simulations for pressures between 1 and 3000 atm are |
| 555 |
|
shown in Fig. \ref{fig:compWater} along with compressibility values |
| 556 |
|
obtained from both other SPC/E simulations and experiment. |
| 529 |
– |
Compressibility values from all references are for applied pressures |
| 530 |
– |
within the range 1 - 1000 atm. |
| 557 |
|
|
| 558 |
|
\begin{figure} |
| 559 |
|
\includegraphics[width=\linewidth]{new_isothermalN} |
| 563 |
|
|
| 564 |
|
Isothermal compressibility values calculated using the number density |
| 565 |
|
(Eq. \ref{eq:BMN}) expression are in good agreement with experimental |
| 566 |
< |
and previous simulation work throughout the 1 - 1000 atm pressure |
| 566 |
> |
and previous simulation work throughout the 1 -- 1000 atm pressure |
| 567 |
|
regime. Compressibilities computed using the Hull volume, however, |
| 568 |
|
deviate dramatically from the experimental values at low applied |
| 569 |
< |
pressures. The reason for this deviation is quite simple; at low |
| 569 |
> |
pressures. The reason for this deviation is quite simple: at low |
| 570 |
|
applied pressures, the liquid is in equilibrium with a vapor phase, |
| 571 |
|
and it is entirely possible for one (or a few) molecules to drift away |
| 572 |
|
from the liquid cluster (see Fig. \ref{fig:coneOfShame}). At low |
| 575 |
|
geometries which include large volumes of empty space. |
| 576 |
|
|
| 577 |
|
\begin{figure} |
| 578 |
< |
\includegraphics[width=\linewidth]{flytest2} |
| 578 |
> |
\includegraphics[width=\linewidth]{coneOfShame} |
| 579 |
|
\caption{At low pressures, the liquid is in equilibrium with the vapor |
| 580 |
|
phase, and isolated molecules can detach from the liquid droplet. |
| 581 |
|
This is expected behavior, but the volume of the convex hull |
| 582 |
< |
includes large regions of empty space. For this reason, |
| 582 |
> |
includes large regions of empty space. For this reason, |
| 583 |
|
compressibilities are computed using local number densities rather |
| 584 |
|
than hull volumes.} |
| 585 |
|
\label{fig:coneOfShame} |
| 596 |
|
different pressures must be done to compute the first derivatives. It |
| 597 |
|
is also possible to compute the compressibility using the fluctuation |
| 598 |
|
dissipation theorem using either fluctuations in the |
| 599 |
< |
volume,\cite{Debenedetti1986}, |
| 599 |
> |
volume,\cite{Debenedetti1986} |
| 600 |
|
\begin{equation} |
| 601 |
|
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle |
| 602 |
|
V \right \rangle ^{2}}{V \, k_{B} \, T}, |
| 603 |
+ |
\label{eq:BMVfluct} |
| 604 |
|
\end{equation} |
| 605 |
|
or, equivalently, fluctuations in the number of molecules within the |
| 606 |
|
fixed region, |
| 607 |
|
\begin{equation} |
| 608 |
|
\kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle |
| 609 |
< |
N \right \rangle ^{2}}{N \, k_{B} \, T}, |
| 609 |
> |
N \right \rangle ^{2}}{N \, k_{B} \, T}. |
| 610 |
> |
\label{eq:BMNfluct} |
| 611 |
|
\end{equation} |
| 612 |
|
Thus, the compressibility of each simulation can be calculated |
| 613 |
< |
entirely independently from all other trajectories. However, the |
| 614 |
< |
resulting compressibilities were still as much as an order of |
| 615 |
< |
magnitude larger than the reference values. However, compressibility |
| 616 |
< |
calculation that relies on the hull volume will suffer these effects. |
| 617 |
< |
WE NEED MORE HERE. |
| 613 |
> |
entirely independently from other trajectories. Compressibility |
| 614 |
> |
calculations that rely on the hull volume will still suffer the |
| 615 |
> |
effects of the empty space due to the vapor phase; for this reason, we |
| 616 |
> |
recommend using the number density (Eq. \ref{eq:BMN}) or number |
| 617 |
> |
density fluctuations (Eq. \ref{eq:BMNfluct}) for computing |
| 618 |
> |
compressibilities. We obtained the results in |
| 619 |
> |
Fig. \ref{fig:compWater} using a sampling radius that was |
| 620 |
> |
approximately 80\% of the mean distance between the center of mass of |
| 621 |
> |
the cluster and the hull atoms. This ratio of sampling radius to |
| 622 |
> |
average hull radius excludes the problematic vapor phase on the |
| 623 |
> |
outside of the cluster while including enough of the liquid phase to |
| 624 |
> |
avoid poor statistics due to fluctuating local densities. |
| 625 |
|
|
| 626 |
+ |
A comparison of the oxygen-oxygen radial distribution functions for |
| 627 |
+ |
SPC/E water simulated using both the Langevin Hull and more |
| 628 |
+ |
traditional periodic boundary methods -- both at 1 atm and 300K -- |
| 629 |
+ |
reveals an understructuring of water in the Langevin Hull that |
| 630 |
+ |
manifests as a slight broadening of the solvation shells. This effect |
| 631 |
+ |
may be due to the introduction of a liquid-vapor interface in the |
| 632 |
+ |
Langevin Hull simulations (an interface which is missing in most |
| 633 |
+ |
periodic simulations of bulk water). Vapor-phase molecules contribute |
| 634 |
+ |
a small but nearly flat portion of the radial distribution function. |
| 635 |
+ |
|
| 636 |
|
\subsection{Molecular orientation distribution at cluster boundary} |
| 637 |
|
|
| 638 |
< |
In order for non-periodic boundary conditions to be widely applicable, |
| 639 |
< |
they must be constructed in such a way that they allow a finite system |
| 640 |
< |
to replicate the properties of the bulk. Early non-periodic |
| 641 |
< |
simulation methods (e.g. hydrophobic boundary potentials) induced |
| 642 |
< |
spurious orientational correlations deep within the simulated |
| 638 |
> |
In order for a non-periodic boundary method to be widely applicable, |
| 639 |
> |
it must be constructed in such a way that they allow a finite system |
| 640 |
> |
to replicate the properties of the bulk. Early non-periodic simulation |
| 641 |
> |
methods (e.g. hydrophobic boundary potentials) induced spurious |
| 642 |
> |
orientational correlations deep within the simulated |
| 643 |
|
system.\cite{Lee1984,Belch1985} This behavior spawned many methods for |
| 644 |
< |
fixing and characterizing the effects of artifical boundaries |
| 644 |
> |
fixing and characterizing the effects of artificial boundaries |
| 645 |
|
including methods which fix the orientations of a set of edge |
| 646 |
|
molecules.\cite{Warshel1978,King1989} |
| 647 |
|
|
| 648 |
|
As described above, the Langevin Hull does not require that the |
| 649 |
|
orientation of molecules be fixed, nor does it utilize an explicitly |
| 650 |
< |
hydrophobic boundary, orientational constraint or radial constraint. |
| 651 |
< |
Therefore, the orientational correlations of the molecules in a water |
| 652 |
< |
cluster are of particular interest in testing this method. Ideally, |
| 653 |
< |
the water molecules on the surface of the cluster will have enough |
| 654 |
< |
mobility into and out of the center of the cluster to maintain a |
| 650 |
> |
hydrophobic boundary, or orientational or radial constraints. |
| 651 |
> |
Therefore, the orientational correlations of the molecules in water |
| 652 |
> |
clusters are of particular interest in testing this method. Ideally, |
| 653 |
> |
the water molecules on the surfaces of the clusters will have enough |
| 654 |
> |
mobility into and out of the center of the cluster to maintain |
| 655 |
|
bulk-like orientational distribution in the absence of orientational |
| 656 |
|
and radial constraints. However, since the number of hydrogen bonding |
| 657 |
|
partners available to molecules on the exterior are limited, it is |
| 658 |
< |
likely that there will be some effective hydrophobicity of the hull. |
| 658 |
> |
likely that there will be an effective hydrophobicity of the hull. |
| 659 |
|
|
| 660 |
< |
To determine the extent of these effects demonstrated by the Langevin |
| 661 |
< |
Hull, we examined the orientationations exhibited by SPC/E water in a |
| 662 |
< |
cluster of 1372 molecules at 300 K and at pressures ranging from 1 - |
| 663 |
< |
1000 atm. The orientational angle of a water molecule is described |
| 660 |
> |
To determine the extent of these effects, we examined the |
| 661 |
> |
orientations exhibited by SPC/E water in a cluster of 1372 |
| 662 |
> |
molecules at 300 K and at pressures ranging from 1 -- 1000 atm. The |
| 663 |
> |
orientational angle of a water molecule is described by |
| 664 |
|
\begin{equation} |
| 665 |
|
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|} |
| 666 |
|
\end{equation} |
| 667 |
|
where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of |
| 668 |
< |
mass and the cluster center of mass and $\vec{\mu}_{i}$ is the vector |
| 669 |
< |
bisecting the H-O-H angle of molecule {\it i} Bulk-like distributions |
| 670 |
< |
will result in $\langle \cos \theta \rangle$ values close to zero. If |
| 671 |
< |
the hull exhibits an overabundance of externally-oriented oxygen sites |
| 672 |
< |
the average orientation will be negative, while dangling hydrogen |
| 673 |
< |
sites will result in positive average orientations. |
| 668 |
> |
mass and the cluster center of mass, and $\vec{\mu}_{i}$ is the vector |
| 669 |
> |
bisecting the H-O-H angle of molecule {\it i}. Bulk-like |
| 670 |
> |
distributions will result in $\langle \cos \theta \rangle$ values |
| 671 |
> |
close to zero. If the hull exhibits an overabundance of |
| 672 |
> |
externally-oriented oxygen sites, the average orientation will be |
| 673 |
> |
negative, while dangling hydrogen sites will result in positive |
| 674 |
> |
average orientations. |
| 675 |
|
|
| 676 |
|
Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values |
| 677 |
|
for molecules in the interior of the cluster (squares) and for |
| 680 |
|
\includegraphics[width=\linewidth]{pAngle} |
| 681 |
|
\caption{Distribution of $\cos{\theta}$ values for molecules on the |
| 682 |
|
interior of the cluster (squares) and for those participating in the |
| 683 |
< |
convex hull (circles) at a variety of pressures. The Langevin hull |
| 683 |
> |
convex hull (circles) at a variety of pressures. The Langevin Hull |
| 684 |
|
exhibits minor dewetting behavior with exposed oxygen sites on the |
| 685 |
|
hull water molecules. The orientational preference for exposed |
| 686 |
|
oxygen appears to be independent of applied pressure. } |
| 692 |
|
orientations. Molecules included in the convex hull show a slight |
| 693 |
|
preference for values of $\cos{\theta} < 0.$ These values correspond |
| 694 |
|
to molecules with oxygen directed toward the exterior of the cluster, |
| 695 |
< |
forming a dangling hydrogen bond acceptor site. |
| 695 |
> |
forming dangling hydrogen bond acceptor sites. |
| 696 |
|
|
| 697 |
< |
In the absence of an electrostatic contribution from the exterior |
| 698 |
< |
bath, the orientational distribution of water molecules included in |
| 699 |
< |
the Langevin Hull will slightly resemble the distribution at a neat |
| 700 |
< |
water liquid/vapor interface. Previous molecular dynamics simulations |
| 701 |
< |
of SPC/E water \cite{Taylor1996} have shown that molecules at the |
| 656 |
< |
liquid/vapor interface favor an orientation where one hydrogen |
| 657 |
< |
protrudes from the liquid phase. This behavior is demonstrated by |
| 658 |
< |
experiments \cite{Du1994} \cite{Scatena2001} showing that |
| 659 |
< |
approximately one-quarter of water molecules at the liquid/vapor |
| 660 |
< |
interface form dangling hydrogen bonds. The negligible preference |
| 661 |
< |
shown in these cluster simulations could be removed through the |
| 662 |
< |
introduction of an implicit solvent model, which would provide the |
| 663 |
< |
missing electrostatic interactions between the cluster molecules and |
| 664 |
< |
the surrounding temperature/pressure bath. |
| 697 |
> |
The orientational preference exhibited by water molecules on the hull |
| 698 |
> |
is significantly weaker than the preference caused by an explicit |
| 699 |
> |
hydrophobic bounding potential. Additionally, the Langevin Hull does |
| 700 |
> |
not require that the orientation of any molecules be fixed in order to |
| 701 |
> |
maintain bulk-like structure, even near the cluster surface. |
| 702 |
|
|
| 703 |
< |
The orientational preference exhibited by hull molecules in the |
| 704 |
< |
Langevin hull is significantly weaker than the preference caused by an |
| 705 |
< |
explicit hydrophobic bounding potential. Additionally, the Langevin |
| 706 |
< |
Hull does not require that the orientation of any molecules be fixed |
| 707 |
< |
in order to maintain bulk-like structure, even at the cluster surface. |
| 703 |
> |
Previous molecular dynamics simulations of SPC/E liquid / vapor |
| 704 |
> |
interfaces using periodic boundary conditions have shown that |
| 705 |
> |
molecules on the liquid side of interface favor a similar orientation |
| 706 |
> |
where oxygen is directed away from the bulk.\cite{Taylor1996} These |
| 707 |
> |
simulations had well-defined liquid and vapor phase regions |
| 708 |
> |
equilibrium and it was observed that {\it vapor} molecules generally |
| 709 |
> |
had one hydrogen protruding from the surface, forming a dangling |
| 710 |
> |
hydrogen bond donor. Our water clusters do not have a true vapor |
| 711 |
> |
region, but rather a few transient molecules that leave the liquid |
| 712 |
> |
droplet (and which return to the droplet relatively quickly). |
| 713 |
> |
Although we cannot obtain an orientational preference of vapor phase |
| 714 |
> |
molecules in a Langevin Hull simulation, but we do agree with previous |
| 715 |
> |
estimates of the orientation of {\it liquid phase} molecules at the |
| 716 |
> |
interface. |
| 717 |
|
|
| 718 |
|
\subsection{Heterogeneous nanoparticle / water mixtures} |
| 719 |
|
|
| 720 |
+ |
To further test the method, we simulated gold nanoparticles ($r = 18$ |
| 721 |
+ |
\AA~, 1433 atoms) solvated by explicit SPC/E water clusters (5000 |
| 722 |
+ |
molecules) using a model for the gold / water interactions that has |
| 723 |
+ |
been used by Dou {\it et. al.} for investigating the separation of |
| 724 |
+ |
water films near hot metal surfaces.\cite{ISI:000167766600035} The |
| 725 |
+ |
Langevin Hull was used to sample pressures of 1, 2, 5, 10, 20, 50, 100 |
| 726 |
+ |
and 200 atm, while all simulations were done at a temperature of 300 |
| 727 |
+ |
K. At these temperatures and pressures, there is no observed |
| 728 |
+ |
separation of the water film from the surface. |
| 729 |
+ |
|
| 730 |
+ |
In Fig. \ref{fig:RhoR} we show the density of water and gold as a |
| 731 |
+ |
function of the distance from the center of the nanoparticle. Higher |
| 732 |
+ |
applied pressures appear to destroy structural correlations in the |
| 733 |
+ |
outermost monolayer of the gold nanoparticle as well as in the water |
| 734 |
+ |
at the near the metal / water interface. Simulations at increased |
| 735 |
+ |
pressures exhibit significant overlap of the gold and water densities, |
| 736 |
+ |
indicating a less well-defined interfacial surface. |
| 737 |
+ |
|
| 738 |
+ |
\begin{figure} |
| 739 |
+ |
\includegraphics[width=\linewidth]{RhoR} |
| 740 |
+ |
\caption{Density profiles of gold and water at the nanoparticle |
| 741 |
+ |
surface. Each curve has been normalized by the average density in |
| 742 |
+ |
the bulk-like region available to the corresponding material. |
| 743 |
+ |
Higher applied pressures de-structure both the gold nanoparticle |
| 744 |
+ |
surface and water at the metal/water interface.} |
| 745 |
+ |
\label{fig:RhoR} |
| 746 |
+ |
\end{figure} |
| 747 |
+ |
|
| 748 |
+ |
At even higher pressures (500 atm and above), problems with the metal |
| 749 |
+ |
- water interaction potential became quite clear. The model we are |
| 750 |
+ |
using appears to have been parameterized for relatively low pressures; |
| 751 |
+ |
it utilizes both shifted Morse and repulsive Morse potentials to model |
| 752 |
+ |
the Au/O and Au/H interactions, respectively. The repulsive wall of |
| 753 |
+ |
the Morse potential does not diverge quickly enough at short distances |
| 754 |
+ |
to prevent water from diffusing into the center of the gold |
| 755 |
+ |
nanoparticles. This behavior is likely not a realistic description of |
| 756 |
+ |
the real physics of the situation. A better model of the gold-water |
| 757 |
+ |
adsorption behavior would require harder repulsive walls to prevent |
| 758 |
+ |
this behavior. |
| 759 |
+ |
|
| 760 |
|
\section{Discussion} |
| 761 |
|
\label{sec:discussion} |
| 762 |
|
|
| 763 |
|
The Langevin Hull samples the isobaric-isothermal ensemble for |
| 764 |
< |
non-periodic systems by coupling the system to an bath characterized |
| 765 |
< |
by pressure, temperature, and solvent viscosity. This enables the |
| 766 |
< |
study of heterogeneous systems composed of materials of significantly |
| 767 |
< |
different compressibilities. Because the boundary is dynamically |
| 768 |
< |
determined during the simulation and the molecules interacting with |
| 769 |
< |
the boundary can change, the method and has minimal perturbations on |
| 770 |
< |
the behavior of molecules at the edges of the simulation. Further |
| 771 |
< |
work on this method will involve implicit electrostatics at the |
| 772 |
< |
boundary (which is missing in the current implementation) as well as |
| 773 |
< |
more sophisticated treatments of the surface geometry (alpha |
| 764 |
> |
non-periodic systems by coupling the system to a bath characterized by |
| 765 |
> |
pressure, temperature, and solvent viscosity. This enables the |
| 766 |
> |
simulation of heterogeneous systems composed of materials with |
| 767 |
> |
significantly different compressibilities. Because the boundary is |
| 768 |
> |
dynamically determined during the simulation and the molecules |
| 769 |
> |
interacting with the boundary can change, the method inflicts minimal |
| 770 |
> |
perturbations on the behavior of molecules at the edges of the |
| 771 |
> |
simulation. Further work on this method will involve implicit |
| 772 |
> |
electrostatics at the boundary (which is missing in the current |
| 773 |
> |
implementation) as well as more sophisticated treatments of the |
| 774 |
> |
surface geometry (alpha |
| 775 |
|
shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight |
| 776 |
|
Cocone\cite{Dey:2003ts}). The non-convex hull geometries are |
| 777 |
|
significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull |
| 782 |
|
|
| 783 |
|
In order to use the Langevin Hull for simulations on parallel |
| 784 |
|
computers, one of the more difficult tasks is to compute the bounding |
| 785 |
< |
surface, facets, and resistance tensors when the processors have |
| 786 |
< |
incomplete information about the entire system's topology. Most |
| 785 |
> |
surface, facets, and resistance tensors when the individual processors |
| 786 |
> |
have incomplete information about the entire system's topology. Most |
| 787 |
|
parallel decomposition methods assign primary responsibility for the |
| 788 |
|
motion of an atomic site to a single processor, and we can exploit |
| 789 |
|
this to efficiently compute the convex hull for the entire system. |
| 821 |
|
The individual hull operations scale with |
| 822 |
|
$\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total |
| 823 |
|
number of sites, and $p$ is the number of processors. These local |
| 824 |
< |
hull operations create a set of $p$ hulls each with approximately |
| 825 |
< |
$\frac{n}{3pr}$ sites (for a cluster of radius $r$). The worst-case |
| 824 |
> |
hull operations create a set of $p$ hulls, each with approximately |
| 825 |
> |
$\frac{n}{3pr}$ sites for a cluster of radius $r$. The worst-case |
| 826 |
|
communication cost for using a ``gather'' operation to distribute this |
| 827 |
|
information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n |
| 828 |
|
\beta (p-1)}{3 r p^2})$, while the final computation of the system |
| 830 |
|
|
| 831 |
|
For a large number of atoms on a moderately parallel machine, the |
| 832 |
|
total costs are dominated by the computations of the individual hulls, |
| 833 |
< |
and communication of these hulls to so the Langevin hull sees roughly |
| 833 |
> |
and communication of these hulls to create the Langevin Hull sees roughly |
| 834 |
|
linear speed-up with increasing processor counts. |
| 835 |
|
|
| 836 |
|
\section*{Acknowledgments} |
| 839 |
|
time was provided by the Center for Research Computing (CRC) at the |
| 840 |
|
University of Notre Dame. |
| 841 |
|
|
| 842 |
+ |
Molecular graphics images were produced using the UCSF Chimera package from |
| 843 |
+ |
the Resource for Biocomputing, Visualization, and Informatics at the |
| 844 |
+ |
University of California, San Francisco (supported by NIH P41 RR001081). |
| 845 |
|
\newpage |
| 846 |
|
|
| 847 |
|
\bibliography{langevinHull} |
