| 18 |
|
\setlength{\belowcaptionskip}{30 pt} |
| 19 |
|
|
| 20 |
|
\bibpunct{[}{]}{,}{s}{}{;} |
| 21 |
< |
\bibliographystyle{aip} |
| 21 |
> |
\bibliographystyle{achemso} |
| 22 |
|
|
| 23 |
|
\begin{document} |
| 24 |
|
|
| 40 |
|
We have developed a new isobaric-isothermal (NPT) algorithm which |
| 41 |
|
applies an external pressure to the facets comprising the convex |
| 42 |
|
hull surrounding the system. A Langevin thermostat is also applied |
| 43 |
< |
to facets of the hull to mimic contact with an external heat |
| 44 |
< |
bath. This new method, the ``Langevin Hull'', performs better than |
| 45 |
< |
traditional affine transform methods for systems containing |
| 46 |
< |
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 |
< |
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. |
| 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. |
| 55 |
|
\end{abstract} |
| 56 |
|
|
| 57 |
|
\newpage |
| 66 |
|
\section{Introduction} |
| 67 |
|
|
| 68 |
|
The most common molecular dynamics methods for sampling configurations |
| 69 |
< |
of an isobaric-isothermal (NPT) ensemble maintain a target pressure in |
| 70 |
< |
a simulation by coupling the volume of the system to a {\it barostat}, |
| 71 |
< |
which is an extra degree of freedom propagated along with the particle |
| 72 |
< |
coordinates. These methods require periodic boundary conditions, |
| 73 |
< |
because when the instantaneous pressure in the system differs from the |
| 74 |
< |
target pressure, the volume is reduced or expanded using {\it affine |
| 75 |
< |
transforms} of the system geometry. An affine transform scales the |
| 76 |
< |
size and shape of the periodic box as well as the particle positions |
| 77 |
< |
within the box (but not the sizes of the particles). The most common |
| 78 |
< |
constant pressure methods, including the Melchionna |
| 79 |
< |
modification\cite{Melchionna1993} to the Nos\'e-Hoover-Andersen |
| 80 |
< |
equations of motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} |
| 81 |
< |
the Berendsen pressure bath,\cite{ISI:A1984TQ73500045} and the |
| 82 |
< |
Langevin Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize |
| 69 |
> |
from an isobaric-isothermal (NPT) ensemble maintain a target pressure |
| 70 |
> |
in a simulation by coupling the volume of the system to a {\it |
| 71 |
> |
barostat}, which is an extra degree of freedom propagated along with |
| 72 |
> |
the particle coordinates. These methods require periodic boundary |
| 73 |
> |
conditions, because when the instantaneous pressure in the system |
| 74 |
> |
differs from the target pressure, the volume is reduced or expanded |
| 75 |
> |
using {\it affine transforms} of the system geometry. An affine |
| 76 |
> |
transform scales the size and shape of the periodic box as well as the |
| 77 |
> |
particle positions within the box (but not the sizes of the |
| 78 |
> |
particles). The most common constant pressure methods, including the |
| 79 |
> |
Melchionna modification\cite{Melchionna1993} to the |
| 80 |
> |
Nos\'e-Hoover-Andersen equations of |
| 81 |
> |
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen |
| 82 |
> |
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin |
| 83 |
> |
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize scaled |
| 84 |
|
coordinate transformation to adjust the box volume. As long as the |
| 85 |
< |
material in the simulation box is essentially a bulk-like liquid which |
| 86 |
< |
has a relatively uniform compressibility, the standard affine |
| 87 |
< |
transform approach provides an excellent way of adjusting the volume |
| 88 |
< |
of the system and applying pressure directly via the interactions |
| 88 |
< |
between atomic sites. |
| 85 |
> |
material in the simulation box has a relatively uniform |
| 86 |
> |
compressibility, the standard affine transform approach provides an |
| 87 |
> |
excellent way of adjusting the volume of the system and applying |
| 88 |
> |
pressure directly via the interactions between atomic sites. |
| 89 |
|
|
| 90 |
|
One problem with this approach appears when the system being simulated |
| 91 |
|
is an inhomogeneous mixture in which portions of the simulation box |
| 100 |
|
|
| 101 |
|
\begin{figure} |
| 102 |
|
\includegraphics[width=\linewidth]{AffineScale2} |
| 103 |
< |
\caption{Affine Scaling constant pressure methods use box-length |
| 104 |
< |
scaling to adjust the volume to adjust to under- or over-pressure |
| 105 |
< |
conditions. In a system with a uniform compressibility (e.g. bulk |
| 106 |
< |
fluids) these methods can work well. In systems containing |
| 107 |
< |
heterogeneous mixtures, the affine scaling moves required to adjust |
| 108 |
< |
the pressure in the high-compressibility regions can cause molecules |
| 109 |
< |
in low compressibility regions to collide.} |
| 103 |
> |
\caption{Affine scaling methods use box-length scaling to adjust the |
| 104 |
> |
volume to adjust to under- or over-pressure conditions. In a system |
| 105 |
> |
with a uniform compressibility (e.g. bulk fluids) these methods can |
| 106 |
> |
work well. In systems containing heterogeneous mixtures, the affine |
| 107 |
> |
scaling moves required to adjust the pressure in the |
| 108 |
> |
high-compressibility regions can cause molecules in low |
| 109 |
> |
compressibility regions to collide.} |
| 110 |
|
\label{affineScale} |
| 111 |
|
\end{figure} |
| 112 |
|
|
| 115 |
|
pressure conditions. The use of periodic boxes to enforce a system |
| 116 |
|
volume requires either 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 |
| 118 |
> |
effect. For example, calculations using typical hydration boxes |
| 119 |
|
solvating a protein under periodic boundary conditions are quite |
| 120 |
< |
expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL |
| 121 |
< |
SIMULATIONS] |
| 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. |
| 128 |
+ |
|
| 129 |
+ |
|
| 130 |
|
\subsection*{Boundary Methods} |
| 131 |
< |
There have been a number of other approaches to explicit |
| 132 |
< |
non-periodicity that focus on constant or nearly-constant {\it volume} |
| 133 |
< |
conditions while maintaining bulk-like behavior. Berkowitz and |
| 134 |
< |
McCammon introduced a stochastic (Langevin) boundary layer inside a |
| 135 |
< |
region of fixed molecules which effectively enforces constant |
| 136 |
< |
temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this |
| 137 |
< |
approach, the stochastic and fixed regions were defined relative to a |
| 138 |
< |
central atom. Brooks and Karplus extended this method to include |
| 139 |
< |
deformable stochastic boundaries.\cite{iii:6312} The stochastic |
| 140 |
< |
boundary approach has been used widely for protein |
| 141 |
< |
simulations. [CITATIONS NEEDED] |
| 131 |
> |
There have been a number of approaches to handle simulations of |
| 132 |
> |
explicitly non-periodic systems that focus on constant or |
| 133 |
> |
nearly-constant {\it volume} conditions while maintaining bulk-like |
| 134 |
> |
behavior. Berkowitz and McCammon introduced a stochastic (Langevin) |
| 135 |
> |
boundary layer inside a region of fixed molecules which effectively |
| 136 |
> |
enforces constant temperature and volume (NVT) |
| 137 |
> |
conditions.\cite{Berkowitz1982} In this approach, the stochastic and |
| 138 |
> |
fixed regions were defined relative to a central atom. Brooks and |
| 139 |
> |
Karplus extended this method to include deformable stochastic |
| 140 |
> |
boundaries.\cite{iii:6312} The stochastic boundary approach has been |
| 141 |
> |
used widely for protein simulations. [CITATIONS NEEDED] |
| 142 |
|
|
| 143 |
|
The electrostatic and dispersive behavior near the boundary has long |
| 144 |
|
been a cause for concern when performing simulations of explicitly |
| 209 |
|
In the following sections, we extend and generalize the approach of |
| 210 |
|
Kohanoff, Caro, and Finnis. The new method, which we are calling the |
| 211 |
|
``Langevin Hull'' applies the external pressure, Langevin drag, and |
| 212 |
< |
random forces on the facets of the {\it hull itself} instead of the |
| 213 |
< |
atomic sites comprising the vertices of the hull. This allows us to |
| 214 |
< |
decouple the external pressure contribution from the drag and random |
| 215 |
< |
force. The methodology is introduced in section \ref{sec:meth}, tests |
| 216 |
< |
on crystalline nanoparticles, liquid clusters, and heterogeneous |
| 217 |
< |
mixtures are detailed in section \ref{sec:tests}. Section |
| 218 |
< |
\ref{sec:discussion} summarizes our findings. |
| 212 |
> |
random forces on the {\it facets of the hull} instead of the atomic |
| 213 |
> |
sites comprising the vertices of the hull. This allows us to decouple |
| 214 |
> |
the external pressure contribution from the drag and random force. |
| 215 |
> |
The methodology is introduced in section \ref{sec:meth}, tests on |
| 216 |
> |
crystalline nanoparticles, liquid clusters, and heterogeneous mixtures |
| 217 |
> |
are detailed in section \ref{sec:tests}. Section \ref{sec:discussion} |
| 218 |
> |
summarizes our findings. |
| 219 |
|
|
| 220 |
|
\section{Methodology} |
| 221 |
|
\label{sec:meth} |
| 245 |
|
simulation. |
| 246 |
|
|
| 247 |
|
\begin{figure} |
| 248 |
< |
\includegraphics[width=\linewidth]{hullSample} |
| 248 |
> |
\includegraphics[width=\linewidth]{solvatedNano} |
| 249 |
|
\caption{The external temperature and pressure bath interacts only |
| 250 |
|
with those atoms on the convex hull (grey surface). The hull is |
| 251 |
< |
computed dynamically at each time step, and molecules dynamically |
| 252 |
< |
move between the interior (Newtonian) region and the Langevin hull.} |
| 251 |
> |
computed dynamically at each time step, and molecules can move |
| 252 |
> |
between the interior (Newtonian) region and the Langevin hull.} |
| 253 |
|
\label{fig:hullSample} |
| 254 |
|
\end{figure} |
| 255 |
|
|
| 370 |
|
\begin{equation} |
| 371 |
|
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
| 372 |
|
\end{equation} |
| 373 |
< |
Note that this treatment explicitly ignores rotations (and |
| 373 |
> |
Note that this treatment ignores rotations (and |
| 374 |
|
translational-rotational coupling) of the facet. In compact systems, |
| 375 |
|
the facets stay relatively fixed in orientation between |
| 376 |
|
configurations, so this appears to be a reasonably good approximation. |
| 380 |
|
molecular dynamics time step, the following process is carried out: |
| 381 |
|
\begin{enumerate} |
| 382 |
|
\item The standard inter-atomic forces ($\nabla_iU$) are computed. |
| 383 |
< |
\item Delaunay triangulation is done using the current atomic |
| 383 |
> |
\item Delaunay triangulation is carried out using the current atomic |
| 384 |
|
configuration. |
| 385 |
|
\item The convex hull is computed and facets are identified. |
| 386 |
|
\item For each facet: |
| 454 |
|
Spohr potential was adopted in depicting the interaction between metal |
| 455 |
|
atoms and the SPC/E water molecules.\cite{ISI:000167766600035} |
| 456 |
|
|
| 457 |
< |
\subsection{Bulk modulus of gold nanoparticles} |
| 457 |
> |
\subsection{Compressibility of gold nanoparticles} |
| 458 |
|
|
| 459 |
< |
The compressibility is well-known for gold, and it provides a good first |
| 460 |
< |
test of how the method compares to other similar methods. |
| 459 |
> |
The compressibility (and its inverse, the bulk modulus) is well-known |
| 460 |
> |
for gold, and is captured well by the embedded atom method |
| 461 |
> |
(EAM)~\cite{PhysRevB.33.7983} potential |
| 462 |
> |
and related multi-body force fields. In particular, the quantum |
| 463 |
> |
Sutton-Chen potential gets nearly quantitative agreement with the |
| 464 |
> |
experimental bulk modulus values, and makes a good first test of how |
| 465 |
> |
the Langevin Hull will perform at large applied pressures. |
| 466 |
|
|
| 467 |
< |
\begin{figure} |
| 468 |
< |
\includegraphics[width=\linewidth]{P_T_combined} |
| 469 |
< |
\caption{Pressure and temperature response of an 18 \AA\ gold |
| 470 |
< |
nanoparticle initially when first placed in the Langevin Hull |
| 471 |
< |
($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa) and starting |
| 467 |
> |
The Sutton-Chen (SC) potentials are based on a model of a metal which |
| 468 |
> |
treats the nuclei and core electrons as pseudo-atoms embedded in the |
| 469 |
> |
electron density due to the valence electrons on all of the other |
| 470 |
> |
atoms in the system.\cite{Chen90} The SC potential has a simple form that closely |
| 471 |
> |
resembles the Lennard Jones potential, |
| 472 |
> |
\begin{equation} |
| 473 |
> |
\label{eq:SCP1} |
| 474 |
> |
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] , |
| 475 |
> |
\end{equation} |
| 476 |
> |
where $V^{pair}_{ij}$ and $\rho_{i}$ are given by |
| 477 |
> |
\begin{equation} |
| 478 |
> |
\label{eq:SCP2} |
| 479 |
> |
V^{pair}_{ij}(r)=\left( \frac{\alpha_{ij}}{r_{ij}}\right)^{n_{ij}}, \rho_{i}=\sum_{j\neq i}\left( \frac{\alpha_{ij}}{r_{ij}}\right) ^{m_{ij}}. |
| 480 |
> |
\end{equation} |
| 481 |
> |
$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for |
| 482 |
> |
interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in |
| 483 |
> |
Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models |
| 484 |
> |
the interactions between the valence electrons and the cores of the |
| 485 |
> |
pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy |
| 486 |
> |
scale, $c_i$ scales the attractive portion of the potential relative |
| 487 |
> |
to the repulsive interaction and $\alpha_{ij}$ is a length parameter |
| 488 |
> |
that assures a dimensionless form for $\rho$. These parameters are |
| 489 |
> |
tuned to various experimental properties such as the density, cohesive |
| 490 |
> |
energy, and elastic moduli for FCC transition metals. The quantum |
| 491 |
> |
Sutton-Chen (QSC) formulation matches these properties while including |
| 492 |
> |
zero-point quantum corrections for different transition |
| 493 |
> |
metals.\cite{PhysRevB.59.3527} |
| 494 |
> |
|
| 495 |
> |
In bulk gold, the experimentally-measured value for the bulk modulus |
| 496 |
> |
is 180.32 GPa, while previous calculations on the QSC potential in |
| 497 |
> |
periodic-boundary simulations of the bulk have yielded values of |
| 498 |
> |
175.53 GPa.\cite{XXX} Using the same force field, we have performed a |
| 499 |
> |
series of relatively short (200 ps) simulations on 40 \r{A} radius |
| 500 |
> |
nanoparticles under the Langevin Hull at a variety of applied |
| 501 |
> |
pressures ranging from 0 GPa to XXX. We obtain a value of 177.547 GPa |
| 502 |
> |
for the bulk modulus for gold using this echnique. |
| 503 |
> |
|
| 504 |
> |
\begin{figure} |
| 505 |
> |
\includegraphics[width=\linewidth]{stacked} |
| 506 |
> |
\caption{The response of the internal pressure and temperature of gold |
| 507 |
> |
nanoparticles when first placed in the Langevin Hull |
| 508 |
> |
($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting |
| 509 |
|
from initial conditions that were far from the bath pressure and |
| 510 |
< |
temperature. The pressure response is rapid, and the thermal |
| 462 |
< |
equilibration depends on both total surface area and the viscosity |
| 463 |
< |
of the bath.} |
| 510 |
> |
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).} |
| 511 |
|
\label{pressureResponse} |
| 512 |
|
\end{figure} |
| 513 |
|
|
| 516 |
|
P}\right) |
| 517 |
|
\end{equation} |
| 518 |
|
|
| 472 |
– |
\begin{figure} |
| 473 |
– |
\includegraphics[width=\linewidth]{compress_tb} |
| 474 |
– |
\caption{Isothermal Compressibility (18 \AA gold nanoparticle)} |
| 475 |
– |
\label{temperatureResponse} |
| 476 |
– |
\end{figure} |
| 477 |
– |
|
| 519 |
|
\subsection{Compressibility of SPC/E water clusters} |
| 520 |
|
|
| 521 |
|
Prior molecular dynamics simulations on SPC/E water (both in |
| 549 |
|
geometries which include large volumes of empty space. |
| 550 |
|
|
| 551 |
|
\begin{figure} |
| 552 |
< |
\includegraphics[width=\linewidth]{flytest2} |
| 552 |
> |
\includegraphics[width=\linewidth]{coneOfShame} |
| 553 |
|
\caption{At low pressures, the liquid is in equilibrium with the vapor |
| 554 |
|
phase, and isolated molecules can detach from the liquid droplet. |
| 555 |
|
This is expected behavior, but the volume of the convex hull |
| 563 |
|
and the hull geometries are much more compact. Because of the |
| 564 |
|
liquid-vapor effect on the convex hull, the regional number density |
| 565 |
|
approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the |
| 566 |
< |
bulk modulus. |
| 566 |
> |
compressibility. |
| 567 |
|
|
| 568 |
|
In both the traditional compressibility formula (Eq. \ref{eq:BM}) and |
| 569 |
|
the number density version (Eq. \ref{eq:BMN}), multiple simulations at |
| 584 |
|
Thus, the compressibility of each simulation can be calculated |
| 585 |
|
entirely independently from all other trajectories. However, the |
| 586 |
|
resulting compressibilities were still as much as an order of |
| 587 |
< |
magnitude larger than the reference values. Any compressibility |
| 587 |
> |
magnitude larger than the reference values. However, compressibility |
| 588 |
|
calculation that relies on the hull volume will suffer these effects. |
| 589 |
|
WE NEED MORE HERE. |
| 590 |
|
|
| 592 |
|
|
| 593 |
|
In order for non-periodic boundary conditions to be widely applicable, |
| 594 |
|
they must be constructed in such a way that they allow a finite system |
| 595 |
< |
to replicate the properties of the bulk. Naturally, this requirement |
| 596 |
< |
has spawned many methods for fixing and characterizing the effects of |
| 597 |
< |
artifical boundaries. Of particular interest regarding the Langevin |
| 598 |
< |
Hull is the orientation of water molecules that are part of the |
| 599 |
< |
geometric hull. Ideally, all molecules in the cluster will have the |
| 600 |
< |
same orientational distribution as bulk water. |
| 595 |
> |
to replicate the properties of the bulk. Early non-periodic |
| 596 |
> |
simulation methods (e.g. hydrophobic boundary potentials) induced |
| 597 |
> |
spurious orientational correlations deep within the simulated |
| 598 |
> |
system.\cite{Lee1984,Belch1985} This behavior spawned many methods for |
| 599 |
> |
fixing and characterizing the effects of artifical boundaries |
| 600 |
> |
including methods which fix the orientations of a set of edge |
| 601 |
> |
molecules.\cite{Warshel1978,King1989} |
| 602 |
|
|
| 603 |
< |
The orientation of molecules at the edges of a simulated cluster has |
| 604 |
< |
long been a concern when performing simulations of explicitly |
| 605 |
< |
non-periodic systems. Early work led to the surface constrained soft |
| 606 |
< |
sphere dipole model (SCSSD) \cite{Warshel1978} in which the surface |
| 607 |
< |
molecules are fixed in a random orientation representative of the bulk |
| 608 |
< |
solvent structural properties. Belch, et al \cite{Belch1985} simulated |
| 609 |
< |
clusters of TIPS2 water surrounded by a hydrophobic bounding |
| 610 |
< |
potential. The spherical hydrophobic boundary induced dangling |
| 611 |
< |
hydrogen bonds at the surface that propagated deep into the cluster, |
| 612 |
< |
affecting 70\% of the 100 molecules in the simulation. This result |
| 613 |
< |
echoes an earlier study which showed that an extended planar |
| 572 |
< |
hydrophobic surface caused orientational preference at the surface |
| 573 |
< |
which extended 7 \r{A} into the liquid simulation cell |
| 574 |
< |
\cite{Lee1984}. The surface constrained all-atom solvent (SCAAS) model |
| 575 |
< |
\cite{King1989} improved upon its SCSSD predecessor. The SCAAS model |
| 576 |
< |
utilizes a polarization constraint which is applied to the surface |
| 577 |
< |
molecules to maintain bulk-like structure at the cluster surface. A |
| 578 |
< |
radial constraint is used to maintain the desired bulk density of the |
| 579 |
< |
liquid. Both constraint forces are applied only to a pre-determined |
| 580 |
< |
number of the outermost molecules. |
| 603 |
> |
As described above, the Langevin Hull does not require that the |
| 604 |
> |
orientation of molecules be fixed, nor does it utilize an explicitly |
| 605 |
> |
hydrophobic boundary, orientational constraint or radial constraint. |
| 606 |
> |
Therefore, the orientational correlations of the molecules in a water |
| 607 |
> |
cluster are of particular interest in testing this method. Ideally, |
| 608 |
> |
the water molecules on the surface of the cluster will have enough |
| 609 |
> |
mobility into and out of the center of the cluster to maintain a |
| 610 |
> |
bulk-like orientational distribution in the absence of orientational |
| 611 |
> |
and radial constraints. However, since the number of hydrogen bonding |
| 612 |
> |
partners available to molecules on the exterior are limited, it is |
| 613 |
> |
likely that there will be some effective hydrophobicity of the hull. |
| 614 |
|
|
| 615 |
< |
In contrast, the Langevin Hull does not require that the orientation |
| 616 |
< |
of molecules be fixed, nor does it utilize an explicitly hydrophobic |
| 617 |
< |
boundary, orientational constraint or radial constraint. The number |
| 618 |
< |
and identity of the molecules included on the convex hull are dynamic |
| 586 |
< |
properties, thus avoiding the formation of an artificial solvent |
| 587 |
< |
boundary layer. The hope is that the water molecules on the surface of |
| 588 |
< |
the cluster, if left to their own devices in the absence of |
| 589 |
< |
orientational and radial constraints, will maintain a bulk-like |
| 590 |
< |
orientational distribution. |
| 591 |
< |
|
| 592 |
< |
To determine the extent of these effects demonstrated by the Langevin Hull, we examined the orientations exhibited by SPC/E water in a cluster of 1372 molecules at 300 K and at pressures ranging from 1 - 1000 atm. |
| 593 |
< |
|
| 594 |
< |
The orientation of a water molecule is described by |
| 595 |
< |
|
| 615 |
> |
To determine the extent of these effects demonstrated by the Langevin |
| 616 |
> |
Hull, we examined the orientationations exhibited by SPC/E water in a |
| 617 |
> |
cluster of 1372 molecules at 300 K and at pressures ranging from 1 - |
| 618 |
> |
1000 atm. The orientational angle of a water molecule is described |
| 619 |
|
\begin{equation} |
| 620 |
|
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|} |
| 621 |
|
\end{equation} |
| 622 |
+ |
where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of |
| 623 |
+ |
mass and the cluster center of mass and $\vec{\mu}_{i}$ is the vector |
| 624 |
+ |
bisecting the H-O-H angle of molecule {\it i} Bulk-like distributions |
| 625 |
+ |
will result in $\langle \cos \theta \rangle$ values close to zero. If |
| 626 |
+ |
the hull exhibits an overabundance of externally-oriented oxygen sites |
| 627 |
+ |
the average orientation will be negative, while dangling hydrogen |
| 628 |
+ |
sites will result in positive average orientations. |
| 629 |
|
|
| 630 |
< |
where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of mass and the cluster center of mass and $\vec{\mu}_{i}$ is the vector bisecting the H-O-H angle of molecule {\it i}. |
| 631 |
< |
|
| 630 |
> |
Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values |
| 631 |
> |
for molecules in the interior of the cluster (squares) and for |
| 632 |
> |
molecules included in the convex hull (circles). |
| 633 |
|
\begin{figure} |
| 603 |
– |
\includegraphics[width=\linewidth]{g_r_theta} |
| 604 |
– |
\caption{Definition of coordinates} |
| 605 |
– |
\label{coords} |
| 606 |
– |
\end{figure} |
| 607 |
– |
|
| 608 |
– |
Fig. 7 shows the probability of each value of $\cos{\theta}$ for molecules in the interior of the cluster (squares) and for molecules included in the convex hull (circles). |
| 609 |
– |
|
| 610 |
– |
\begin{figure} |
| 634 |
|
\includegraphics[width=\linewidth]{pAngle} |
| 635 |
< |
\caption{SPC/E water clusters: only minor dewetting at the boundary} |
| 636 |
< |
\label{pAngle} |
| 635 |
> |
\caption{Distribution of $\cos{\theta}$ values for molecules on the |
| 636 |
> |
interior of the cluster (squares) and for those participating in the |
| 637 |
> |
convex hull (circles) at a variety of pressures. The Langevin hull |
| 638 |
> |
exhibits minor dewetting behavior with exposed oxygen sites on the |
| 639 |
> |
hull water molecules. The orientational preference for exposed |
| 640 |
> |
oxygen appears to be independent of applied pressure. } |
| 641 |
> |
\label{fig:pAngle} |
| 642 |
|
\end{figure} |
| 643 |
|
|
| 644 |
< |
As expected, interior molecules (those not included in the convex hull) maintain a bulk-like structure with a uniform distribution of orientations. Molecules included in the convex hull show a slight preference for values of $\cos{\theta} < 0.$ These values correspond to molecules with a hydrogen directed toward the exterior of the cluster, forming a dangling hydrogen bond. |
| 644 |
> |
As expected, interior molecules (those not included in the convex |
| 645 |
> |
hull) maintain a bulk-like structure with a uniform distribution of |
| 646 |
> |
orientations. Molecules included in the convex hull show a slight |
| 647 |
> |
preference for values of $\cos{\theta} < 0.$ These values correspond |
| 648 |
> |
to molecules with oxygen directed toward the exterior of the cluster, |
| 649 |
> |
forming a dangling hydrogen bond acceptor site. |
| 650 |
|
|
| 651 |
< |
In the absence of an electrostatic contribution from the exterior bath, the orientational distribution of water molecules included in the Langevin Hull will slightly resemble the distribution at a neat water liquid/vapor interface. Previous molecular dynamics simulations of SPC/E water \cite{Taylor1996} have shown that molecules at the liquid/vapor interface favor an orientation where one hydrogen protrudes from the liquid phase. This behavior is demonstrated by experiments \cite{Du1994} \cite{Scatena2001} showing that approximately one-quarter of water molecules at the liquid/vapor interface form dangling hydrogen bonds. The negligible preference shown in these cluster simulations could be removed through the introduction of an implicit solvent model, which would provide the missing electrostatic interactions between the cluster molecules and the surrounding temperature/pressure bath. |
| 651 |
> |
In the absence of an electrostatic contribution from the exterior |
| 652 |
> |
bath, the orientational distribution of water molecules included in |
| 653 |
> |
the Langevin Hull will slightly resemble the distribution at a neat |
| 654 |
> |
water liquid/vapor interface. Previous molecular dynamics simulations |
| 655 |
> |
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. |
| 665 |
|
|
| 666 |
< |
The orientational preference exhibited by hull molecules is significantly weaker than the preference caused by an explicit hydrophobic bounding potential. Additionally, the Langevin Hull does not require that the orientation of any molecules be fixed in order to maintain bulk-like structure, even at the cluster surface. |
| 666 |
> |
The orientational preference exhibited by hull molecules in the |
| 667 |
> |
Langevin hull is significantly weaker than the preference caused by an |
| 668 |
> |
explicit hydrophobic bounding potential. Additionally, the Langevin |
| 669 |
> |
Hull does not require that the orientation of any molecules be fixed |
| 670 |
> |
in order to maintain bulk-like structure, even at the cluster surface. |
| 671 |
|
|
| 672 |
|
\subsection{Heterogeneous nanoparticle / water mixtures} |
| 673 |
|
|
| 674 |
|
\section{Discussion} |
| 675 |
|
\label{sec:discussion} |
| 676 |
|
|
| 677 |
+ |
The Langevin Hull samples the isobaric-isothermal ensemble for |
| 678 |
+ |
non-periodic systems by coupling the system to an bath characterized |
| 679 |
+ |
by pressure, temperature, and solvent viscosity. This enables the |
| 680 |
+ |
study of heterogeneous systems composed of materials of significantly |
| 681 |
+ |
different compressibilities. Because the boundary is dynamically |
| 682 |
+ |
determined during the simulation and the molecules interacting with |
| 683 |
+ |
the boundary can change, the method and has minimal perturbations on |
| 684 |
+ |
the behavior of molecules at the edges of the simulation. Further |
| 685 |
+ |
work on this method will involve implicit electrostatics at the |
| 686 |
+ |
boundary (which is missing in the current implementation) as well as |
| 687 |
+ |
more sophisticated treatments of the surface geometry (alpha |
| 688 |
+ |
shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight |
| 689 |
+ |
Cocone\cite{Dey:2003ts}). The non-convex hull geometries are |
| 690 |
+ |
significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull |
| 691 |
+ |
($\mathcal{O}(N \log N)$), but would enable the use of hull volumes |
| 692 |
+ |
directly in computing the compressibility of the sample. |
| 693 |
+ |
|
| 694 |
|
\section*{Appendix A: Computing Convex Hulls on Parallel Computers} |
| 695 |
|
|
| 696 |
+ |
In order to use the Langevin Hull for simulations on parallel |
| 697 |
+ |
computers, one of the more difficult tasks is to compute the bounding |
| 698 |
+ |
surface, facets, and resistance tensors when the processors have |
| 699 |
+ |
incomplete information about the entire system's topology. Most |
| 700 |
+ |
parallel decomposition methods assign primary responsibility for the |
| 701 |
+ |
motion of an atomic site to a single processor, and we can exploit |
| 702 |
+ |
this to efficiently compute the convex hull for the entire system. |
| 703 |
+ |
|
| 704 |
+ |
The basic idea involves splitting the point cloud into |
| 705 |
+ |
spatially-overlapping subsets and computing the convex hulls for each |
| 706 |
+ |
of the subsets. The points on the convex hull of the entire system |
| 707 |
+ |
are all present on at least one of the subset hulls. The algorithm |
| 708 |
+ |
works as follows: |
| 709 |
+ |
\begin{enumerate} |
| 710 |
+ |
\item Each processor computes the convex hull for its own atomic sites |
| 711 |
+ |
(left panel in Fig. \ref{fig:parallel}). |
| 712 |
+ |
\item The Hull vertices from each processor are communicated to all of |
| 713 |
+ |
the processors, and each processor assembles a complete list of hull |
| 714 |
+ |
sites (this is much smaller than the original number of points in |
| 715 |
+ |
the point cloud). |
| 716 |
+ |
\item Each processor computes the global convex hull (right panel in |
| 717 |
+ |
Fig. \ref{fig:parallel}) using only those points that are the union |
| 718 |
+ |
of sites gathered from all of the subset hulls. Delaunay |
| 719 |
+ |
triangulation is then done to obtain the facets of the global hull. |
| 720 |
+ |
\end{enumerate} |
| 721 |
+ |
|
| 722 |
+ |
\begin{figure} |
| 723 |
+ |
\includegraphics[width=\linewidth]{parallel} |
| 724 |
+ |
\caption{When the sites are distributed among many nodes for parallel |
| 725 |
+ |
computation, the processors first compute the convex hulls for their |
| 726 |
+ |
own sites (dashed lines in left panel). The positions of the sites |
| 727 |
+ |
that make up the subset hulls are then communicated to all |
| 728 |
+ |
processors (middle panel). The convex hull of the system (solid line in |
| 729 |
+ |
right panel) is the convex hull of the points on the union of the subset |
| 730 |
+ |
hulls.} |
| 731 |
+ |
\label{fig:parallel} |
| 732 |
+ |
\end{figure} |
| 733 |
+ |
|
| 734 |
+ |
The individual hull operations scale with |
| 735 |
+ |
$\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total |
| 736 |
+ |
number of sites, and $p$ is the number of processors. These local |
| 737 |
+ |
hull operations create a set of $p$ hulls each with approximately |
| 738 |
+ |
$\frac{n}{3pr}$ sites (for a cluster of radius $r$). The worst-case |
| 739 |
+ |
communication cost for using a ``gather'' operation to distribute this |
| 740 |
+ |
information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n |
| 741 |
+ |
\beta (p-1)}{3 r p^2})$, while the final computation of the system |
| 742 |
+ |
hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$. |
| 743 |
+ |
|
| 744 |
+ |
For a large number of atoms on a moderately parallel machine, the |
| 745 |
+ |
total costs are dominated by the computations of the individual hulls, |
| 746 |
+ |
and communication of these hulls to so the Langevin hull sees roughly |
| 747 |
+ |
linear speed-up with increasing processor counts. |
| 748 |
+ |
|
| 749 |
|
\section*{Acknowledgments} |
| 750 |
|
Support for this project was provided by the |
| 751 |
|
National Science Foundation under grant CHE-0848243. Computational |
| 752 |
|
time was provided by the Center for Research Computing (CRC) at the |
| 753 |
|
University of Notre Dame. |
| 754 |
|
|
| 755 |
+ |
Molecular graphics images were produced using the UCSF Chimera package from |
| 756 |
+ |
the Resource for Biocomputing, Visualization, and Informatics at the |
| 757 |
+ |
University of California, San Francisco (supported by NIH P41 RR001081). |
| 758 |
|
\newpage |
| 759 |
|
|
| 760 |
|
\bibliography{langevinHull} |
