trunk/langevinHull/langevinHull.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
\usepackage{graphicx}
\usepackage{multirow}
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures
%\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

\bibpunct{}{}{,}{s}{}{;}
\bibliographystyle{achemso}

\begin{document}

\title{The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems}

\author{Charles F. Vardeman II, Kelsey M. Stocker, and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}
  We have developed a new isobaric-isothermal (NPT) algorithm which
  applies an external pressure to the facets comprising the convex
  hull surrounding the system.  A Langevin thermostat is also applied
  to the facets to mimic contact with an external heat bath. This new
  method, the ``Langevin Hull'', can handle heterogeneous mixtures of
  materials with different compressibilities.  These systems are
  problematic for traditional affine transform methods.  The Langevin
  Hull does not suffer from the edge effects of boundary potential
  methods, and allows realistic treatment of both external pressure
  and thermal conductivity due to the presence of an implicit solvent.
  We apply this method to several different systems including bare
  metal nanoparticles, nanoparticles in an explicit solvent, as well
  as clusters of liquid water. The predicted mechanical properties of
  these systems are in good agreement with experimental data and
  previous simulation work.
\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Introduction}

The most common molecular dynamics methods for sampling configurations
from an isobaric-isothermal (NPT) ensemble maintain a target pressure
in a simulation by coupling the volume of the system to a {\it
  barostat}, which is an extra degree of freedom propagated along with
the particle coordinates.  These methods require periodic boundary
conditions, because when the instantaneous pressure in the system
differs from the target pressure, the volume is reduced or expanded
using {\it affine transforms} of the system geometry. An affine
transform scales the size and shape of the periodic box as well as the
particle positions within the box (but not the sizes of the
particles). The most common constant pressure methods, including the
Melchionna modification\cite{Melchionna1993} to the
Nos\'e-Hoover-Andersen equations of
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize scaled
coordinate transformation to adjust the box volume.  As long as the
material in the simulation box has a relatively uniform
compressibility, the standard affine transform approach provides an
excellent way of adjusting the volume of the system and applying
pressure directly via the interactions between atomic sites.

One problem with this approach appears when the system being simulated
is an inhomogeneous mixture in which portions of the simulation box
are incompressible relative to other portions.  Examples include
simulations of metallic nanoparticles in liquid environments, proteins
at ice / water interfaces, as well as other heterogeneous or
interfacial environments.  In these cases, the affine transform of
atomic coordinates will either cause numerical instability when the
sites in the incompressible medium collide with each other, or will
lead to inefficient sampling of system volumes if the barostat is set
slow enough to avoid the instabilities in the incompressible region.

\begin{figure}
\includegraphics[width=\linewidth]{AffineScale2}
\caption{Affine scaling methods use box-length scaling to adjust the
  volume to adjust to under- or over-pressure conditions. In a system
  with a uniform compressibility (e.g. bulk fluids) these methods can
  work well.  In systems containing heterogeneous mixtures, the affine
  scaling moves required to adjust the pressure in the
  high-compressibility regions can cause molecules in low
  compressibility regions to collide.}
\label{affineScale}
\end{figure}

One may also wish to avoid affine transform periodic boundary methods
to simulate {\it explicitly non-periodic systems} under constant
pressure conditions. The use of periodic boxes to enforce a system
volume requires either effective solute concentrations that are much
higher than desirable, or unreasonable system sizes to avoid this
effect.  For example, calculations using typical hydration boxes
solvating a protein under periodic boundary conditions are quite
expensive.  A 62 \AA$^3$ box of water solvating a moderately small
protein like hen egg white lysozyme (PDB code: 1LYZ) yields an
effective protein concentration of 100 mg/mL.\cite{Asthagiri20053300}

{\it Total} protein concentrations in the cell are typically on the
order of 160-310 mg/ml,\cite{Brown1991195} and individual proteins
have concentrations orders of magnitude lower than this in the
cellular environment. The effective concentrations of single proteins
in simulations may have significant effects on the structure and
dynamics of simulated systems.

\subsection*{Boundary Methods}
There have been a number of approaches to handle simulations of
explicitly non-periodic systems that focus on constant or
nearly-constant {\it volume} conditions while maintaining bulk-like
behavior.  Berkowitz and McCammon introduced a stochastic (Langevin)
boundary layer inside a region of fixed molecules which effectively
enforces constant temperature and volume (NVT)
conditions.\cite{Berkowitz1982} In this approach, the stochastic and
fixed regions were defined relative to a central atom.  Brooks and
Karplus extended this method to include deformable stochastic
boundaries.\cite{iii:6312} The stochastic boundary approach has been
used widely for protein simulations.

The electrostatic and dispersive behavior near the boundary has long
been a cause for concern when performing simulations of explicitly
non-periodic systems.  Early work led to the surface constrained soft
sphere dipole model (SCSSD)\cite{Warshel1978} in which the surface
molecules are fixed in a random orientation representative of the bulk
solvent structural properties. Belch {\it et al.}\cite{Belch1985}
simulated clusters of TIPS2 water surrounded by a hydrophobic bounding
potential. The spherical hydrophobic boundary induced dangling
hydrogen bonds at the surface that propagated deep into the cluster,
affecting most of the molecules in the simulation.  This result echoes
an earlier study which showed that an extended planar hydrophobic
surface caused orientational preferences at the surface which extended
relatively deep (7 \AA) into the liquid simulation cell.\cite{Lee1984}
The surface constrained all-atom solvent (SCAAS) model \cite{King1989}
improved upon its SCSSD predecessor. The SCAAS model utilizes a
polarization constraint which is applied to the surface molecules to
maintain bulk-like structure at the cluster surface. A radial
constraint is used to maintain the desired bulk density of the
liquid. Both constraint forces are applied only to a pre-determined
number of the outermost molecules.

Beglov and Roux have developed a boundary model in which the hard
sphere boundary has a radius that varies with the instantaneous
configuration of the solute (and solvent) molecules.\cite{beglov:9050}
This model contains a clear pressure and surface tension contribution
to the free energy.

\subsection*{Restraining Potentials}
Restraining {\it potentials} introduce repulsive potentials at the
surface of a sphere or other geometry.  The solute and any explicit
solvent are therefore restrained inside the range defined by the
external potential.  Often the potentials include a weak short-range
attraction to maintain the correct density at the boundary.  Beglov
and Roux have also introduced a restraining boundary potential which
relaxes dynamically depending on the solute geometry and the force the
explicit system exerts on the shell.\cite{Beglov:1995fk}

Recently, Krilov {\it et al.} introduced a {\it flexible} boundary
model that uses a Lennard-Jones potential between the solvent
molecules and a boundary which is determined dynamically from the
position of the nearest solute atom.\cite{LiY._jp046852t,Zhu:2008fk} This
approach allows the confining potential to prevent solvent molecules
from migrating too far from the solute surface, while providing a weak
attractive force pulling the solvent molecules towards a fictitious
bulk solvent.  Although this approach is appealing and has physical
motivation, nanoparticles do not deform far from their original
geometries even at temperatures which vaporize the nearby solvent. For
the systems like this, the flexible boundary model will be nearly
identical to a fixed-volume restraining potential.

\subsection*{Hull methods}
The approach of Kohanoff, Caro, and Finnis is the most promising of
the methods for introducing both constant pressure and temperature
into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru}
This method is based on standard Langevin dynamics, but the Brownian
or random forces are allowed to act only on peripheral atoms and exert
forces in a direction that is inward-facing relative to the facets of
a closed bounding surface.  The statistical distribution of the random
forces are uniquely tied to the pressure in the external reservoir, so
the method can be shown to sample the isobaric-isothermal ensemble.
Kohanoff {\it et al.} used a Delaunay tessellation to generate a
bounding surface surrounding the outermost atoms in the simulated
system.  This is not the only possible triangulated outer surface, but
guarantees that all of the random forces point inward towards the
cluster.

In the following sections, we extend and generalize the approach of
Kohanoff, Caro, and Finnis. The new method, which we are calling the
``Langevin Hull'' applies the external pressure, Langevin drag, and
random forces on the {\it facets of the hull} instead of the atomic
sites comprising the vertices of the hull.  This allows us to decouple
the external pressure contribution from the drag and random force.
The methodology is introduced in section \ref{sec:meth}, tests on
crystalline nanoparticles, liquid clusters, and heterogeneous mixtures
are detailed in section \ref{sec:tests}.  Section \ref{sec:discussion}
summarizes our findings.

\section{Methodology}
\label{sec:meth}

The Langevin Hull uses an external bath at a fixed constant pressure
($P$) and temperature ($T$) with an effective solvent viscosity
($\eta$).  This bath interacts only with the objects on the exterior
hull of the system.  Defining the hull of the atoms in a simulation is
done in a manner similar to the approach of Kohanoff, Caro and
Finnis.\cite{Kohanoff:2005qm} That is, any instantaneous configuration
of the atoms in the system is considered as a point cloud in three
dimensional space.  Delaunay triangulation is used to find all facets
between coplanar
neighbors.\cite{delaunay,springerlink:10.1007/BF00977785} In highly
symmetric point clouds, facets can contain many atoms, but in all but
the most symmetric of cases, the facets are simple triangles in
3-space which contain exactly three atoms.

The convex hull is the set of facets that have {\it no concave
  corners} at an atomic site.\cite{Barber96,EDELSBRUNNER:1994oq} This
eliminates all facets on the interior of the point cloud, leaving only
those exposed to the bath. Sites on the convex hull are dynamic; as
molecules re-enter the cluster, all interactions between atoms on that
molecule and the external bath are removed.  Since the edge is
determined dynamically as the simulation progresses, no {\it a priori}
geometry is defined. The pressure and temperature bath interacts only
with the atoms on the edge and not with atoms interior to the
simulation.

\begin{figure}
\includegraphics[width=\linewidth]{solvatedNano}
\caption{The external temperature and pressure bath interacts only
  with those atoms on the convex hull (grey surface).  The hull is
  computed dynamically at each time step, and molecules can move 
  between the interior (Newtonian) region and the Langevin Hull.}
\label{fig:hullSample}
\end{figure}

Atomic sites in the interior of the simulation move under standard
Newtonian dynamics,
\begin{equation}
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U,
\label{eq:Newton}
\end{equation}
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the
instantaneous velocity of site $i$ at time $t$, and $U$ is the total
potential energy.  For atoms on the exterior of the cluster
(i.e. those that occupy one of the vertices of the convex hull), the
equation of motion is modified with an external force, ${\mathbf
  F}_i^{\mathrm ext}$:
\begin{equation}
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}.
\end{equation}

The external bath interacts indirectly with the atomic sites through
the intermediary of the hull facets.  Since each vertex (or atom)
provides one corner of a triangular facet, the force on the facets are
divided equally to each vertex.  However, each vertex can participate
in multiple facets, so the resultant force is a sum over all facets
$f$ containing vertex $i$:
\begin{equation}
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\
    } f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\  {\mathbf
  F}_f^{\mathrm ext}
\end{equation}

The external pressure bath applies a force to the facets of the convex
hull in direct proportion to the area of the facet, while the thermal
coupling depends on the solvent temperature, viscosity and the size
and shape of each facet. The thermal interactions are expressed as a
standard Langevin description of the forces,
\begin{equation}
\begin{array}{rclclcl}
{\mathbf F}_f^{\text{ext}} & = &  \text{external pressure} & + & \text{drag force} & + & \text{random force} \\
& = &  -\hat{n}_f P A_f  & - & \Xi_f(t) {\mathbf v}_f(t)  & + & {\mathbf R}_f(t)
\end{array}
\end{equation}
Here, $A_f$ and $\hat{n}_f$ are the area and (outward-facing) normal
vectors for facet $f$, respectively.  ${\mathbf v}_f(t)$ is the
velocity of the facet centroid,
\begin{equation}
{\mathbf v}_f(t) =  \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i,
\end{equation}
and $\Xi_f(t)$ is an approximate ($3 \times 3$) resistance tensor that
depends on the geometry and surface area of facet $f$ and the
viscosity of the bath.  The resistance tensor is related to the
fluctuations of the random force, $\mathbf{R}(t)$, by the
fluctuation-dissipation theorem,
\begin{eqnarray}
\left< {\mathbf R}_f(t) \right> & = & 0 \\
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\
\Xi_f(t)\delta(t-t^\prime).
\label{eq:randomForce}
\end{eqnarray}

Once the resistance tensor is known for a given facet, a stochastic
vector that has the properties in Eq. (\ref{eq:randomForce}) can be
calculated efficiently by carrying out a Cholesky decomposition to
obtain the square root matrix of the resistance tensor,
\begin{equation} 
\Xi_f = {\bf S} {\bf S}^{T},
\label{eq:Cholesky}
\end{equation} 
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
vector with the statistics required for the random force can then be
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which
has elements chosen from a Gaussian distribution, such that:
\begin{equation}
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
\end{equation}
where $\delta t$ is the timestep in use during the simulation. The
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to
have the correct properties required by Eq. (\ref{eq:randomForce}).

Our treatment of the resistance tensor is approximate.  $\Xi_f$ for a
rigid triangular plate would normally be treated as a $6 \times 6$
tensor that includes translational and rotational drag as well as
translational-rotational coupling. The computation of resistance
tensors for rigid bodies has been detailed
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun:2008fk}
but the standard approach involving bead approximations would be
prohibitively expensive if it were recomputed at each step in a
molecular dynamics simulation.

Instead, we are utilizing an approximate resistance tensor obtained by
first constructing the Oseen tensor for the interaction of the
centroid of the facet ($f$) with each of the subfacets $\ell=1,2,3$,
\begin{equation}
T_{\ell f}=\frac{A_\ell}{8\pi\eta R_{\ell f}}\left(I +
  \frac{\mathbf{R}_{\ell f}\mathbf{R}_{\ell f}^T}{R_{\ell f}^2}\right)
\end{equation}
Here, $A_\ell$ is the area of subfacet $\ell$ which is a triangle
containing two of the vertices of the facet along with the centroid.
$\mathbf{R}_{\ell f}$ is the vector between the centroid of facet $f$
and the centroid of sub-facet $\ell$, and $I$ is the ($3 \times 3$)
identity matrix.  $\eta$ is the viscosity of the external bath.

\begin{figure}
\includegraphics[width=\linewidth]{hydro}
\caption{The resistance tensor $\Xi$ for a facet comprising sites $i$,
  $j$, and $k$ is constructed using Oseen tensor contributions between
  the centoid of the facet $f$ and each of the sub-facets ($i,f,j$),
  ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets are
  located at $1$, $2$, and $3$, and the area of each sub-facet is
  easily computed using half the cross product of two of the edges.}
\label{hydro}
\end{figure}

The tensors for each of the sub-facets are added together, and the
resulting matrix is inverted to give a $3 \times 3$ resistance tensor
for translations of the triangular facet,
\begin{equation}
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}.
\end{equation}
Note that this treatment ignores rotations (and
translational-rotational coupling) of the facet.  In compact systems,
the facets stay relatively fixed in orientation between
configurations, so this appears to be a reasonably good approximation.

We have implemented this method by extending the Langevin dynamics
integrator in our code, OpenMD.\cite{Meineke2005,open_md}  At each
molecular dynamics time step, the following process is carried out:
\begin{enumerate}
\item The standard inter-atomic forces ($\nabla_iU$) are computed.
\item Delaunay triangulation is carried out using the current atomic
  configuration.
\item The convex hull is computed and facets are identified.
\item For each facet:
\begin{itemize}
\item[a.] The force from the pressure bath ($-\hat{n}_fPA_f$) is
  computed.
\item[b.] The resistance tensor ($\Xi_f(t)$) is computed using the
  viscosity ($\eta$) of the bath.
\item[c.] Facet drag ($-\Xi_f(t) \mathbf{v}_f(t)$) forces are
  computed.
\item[d.] Random forces ($\mathbf{R}_f(t)$) are computed using the
  resistance tensor and the temperature ($T$) of the bath.
\end{itemize}
\item The facet forces are divided equally among the vertex atoms.
\item Atomic positions and velocities are propagated.
\end{enumerate}
The Delaunay triangulation and computation of the convex hull are done
using calls to the qhull library.\cite{Q_hull} There is a minimal
penalty for computing the convex hull and resistance tensors at each
step in the molecular dynamics simulation (roughly 0.02 $\times$ cost
of a single force evaluation), and the convex hull is remarkably easy
to parallelize on distributed memory machines (see Appendix A).

\section{Tests \& Applications}
\label{sec:tests}

To test the new method, we have carried out simulations using the
Langevin Hull on: 1) a crystalline system (gold nanoparticles), 2) a
liquid droplet (SPC/E water),\cite{Berendsen1987} and 3) a
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,
\begin{equation}
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right
)_{T}.
\label{eq:BM}
\end{equation}

One problem with eliminating periodic boundary conditions and
simulation boxes is that the volume of a three-dimensional point cloud
is not well-defined. In order to compute the compressibility of a
bulk material, we make an assumption that the number density, $\rho =
\frac{N}{V}$, is uniform within some region of the point cloud. The
compressibility can then be expressed in terms of the average number
of particles in that region,
\begin{equation}
\kappa_{T} = -\frac{1}{N} \left ( \frac{\partial N}{\partial P} \right
)_{T}.
\label{eq:BMN}
\end{equation}
The region we used is a spherical volume of 20 \AA\ radius centered in
the middle of the cluster with a roughly 25 \AA\ radius. $N$ is the average number of molecules
found within this region throughout a given simulation. The geometry
of the region is arbitrary, and any bulk-like portion of the
cluster can be used to compute the compressibility. 

One might assume that the volume of the convex hull could simply be
taken as the system volume $V$ in the compressibility expression
(Eq. \ref{eq:BM}), but this has implications at lower pressures (which
are explored in detail in the section on water droplets).

The metallic force field in use for the gold nanoparticles is the
quantum Sutton-Chen (QSC) model.\cite{PhysRevB.59.3527} In all
simulations involving point charges, we utilized damped shifted-force
(DSF) electrostatics\cite{Fennell06} which is a variant of the Wolf
summation\cite{wolf:8254} that has been shown to provide good forces
and torques on molecular models for water in a computationally
efficient manner.\cite{Fennell06} The damping parameter ($\alpha$) was
set to 0.18 \AA$^{-1}$, and the cutoff radius was set to 12 \AA.  The
Spohr potential was adopted in depicting the interaction between metal
atoms and the SPC/E water molecules.\cite{ISI:000167766600035}

\subsection{Bulk Modulus of gold nanoparticles}

The compressibility (and its inverse, the bulk modulus) is well-known
for gold, and is captured well by the embedded atom method
(EAM)~\cite{PhysRevB.33.7983} potential and related multi-body force
fields.  In particular, the quantum Sutton-Chen potential gets nearly
quantitative agreement with the experimental bulk modulus values, and
makes a good first test of how the Langevin Hull will perform at large
applied pressures.

The Sutton-Chen (SC) potentials are based on a model of a metal which
treats the nuclei and core electrons as pseudo-atoms embedded in the
electron density due to the valence electrons on all of the other
atoms in the system.\cite{Chen90} The SC potential has a simple form
that closely resembles the Lennard Jones potential,
\begin{equation}
\label{eq:SCP1} 
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] , 
\end{equation}
where $V^{pair}_{ij}$ and $\rho_{i}$ are given by 
\begin{equation}
\label{eq:SCP2} 
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}}. 
\end{equation}
$V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for
interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in
Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models
the interactions between the valence electrons and the cores of the
pseudo-atoms. $D_{ij}$ and $D_{ii}$ set the appropriate overall energy
scale, $c_i$ scales the attractive portion of the potential relative
to the repulsive interaction and $\alpha_{ij}$ is a length parameter
that assures a dimensionless form for $\rho$. These parameters are
tuned to various experimental properties such as the density, cohesive
energy, and elastic moduli for FCC transition metals. The quantum
Sutton-Chen (QSC) formulation matches these properties while including
zero-point quantum corrections for different transition
metals.\cite{PhysRevB.59.3527,QSC2}

In bulk gold, the experimentally-measured value for the bulk modulus
is 180.32 GPa, while previous calculations on the QSC potential in
periodic-boundary simulations of the bulk crystal have yielded values
of 175.53 GPa.\cite{QSC2} Using the same force field, we have
performed a series of 1 ns simulations on gold nanoparticles of three
different radii: 20 \AA~ (1985 atoms), 30 \AA~ (6699 atoms), and 40
\AA~ (15707 atoms) utilizing the Langevin Hull at a variety of applied
pressures ranging from 0 -- 10 GPa.  For the 40 \AA~ radius
nanoparticle we obtain a value of 177.55 GPa for the bulk modulus of
gold, in close agreement with both previous simulations and the
experimental bulk modulus reported for gold single
crystals.\cite{Collard1991} The smaller gold nanoparticles (30 and 20
\AA~ radii) have calculated bulk moduli of 215.58 and 208.86 GPa,
respectively, indicating that smaller nanoparticles are somewhat
stiffer (less compressible) than the larger nanoparticles.  This
stiffening of the small nanoparticles may be related to their high
degree of surface curvature, resulting in a lower coordination number
of surface atoms relative to the the surface atoms in the 40 \AA~
radius particle.

We obtain a gold lattice constant of 4.051 \AA~ using the Langevin
Hull at 1 atm, close to the experimentally-determined value for bulk
gold and the value for gold simulated using the QSC potential and
periodic boundary conditions (4.079 \AA~ and 4.088\AA~,
respectively).\cite{QSC2} The slightly smaller calculated lattice
constant is most likely due to the presence of surface tension in the
non-periodic Langevin Hull cluster, an effect absent from a bulk
simulation. The specific heat of a 40 \AA~ gold nanoparticle under the
Langevin Hull at 1 atm is 24.914 $\mathrm {\frac{J}{mol \, K}}$, which
compares very well with the experimental value of 25.42 $\mathrm
{\frac{J}{mol \, K}}$.

\begin{figure}
\includegraphics[width=\linewidth]{stacked}
\caption{The response of the internal pressure and temperature of gold
  nanoparticles when first placed in the Langevin Hull
  ($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting
  from initial conditions that were far from the bath pressure and
  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).}
\label{fig:pressureResponse}
\end{figure}

We note that the Langevin Hull produces rapidly-converging behavior
for structures that are started far from equilibrium.  In
Fig. \ref{fig:pressureResponse} we show how the pressure and
temperature respond to the Langevin Hull for nanoparticles that were
initialized far from the target pressure and temperature.  As
expected, the rate at which thermal equilibrium is achieved depends on
the total surface area of the cluster exposed to the bath as well as
the bath viscosity.  Pressure that is applied suddenly to a cluster
can excite breathing vibrations, but these rapidly damp out (on time
scales of 30 -- 50 ps).

\subsection{Compressibility of SPC/E water clusters}

Prior molecular dynamics simulations on SPC/E water (both in
NVT~\cite{Glattli2002} and NPT~\cite{Motakabbir1990, Pi2009}
ensembles) have yielded values for the isothermal compressibility that
agree well with experiment.\cite{Fine1973} The results of two
different approaches for computing the isothermal compressibility from
Langevin Hull simulations for pressures between 1 and 3000 atm are
shown in Fig. \ref{fig:compWater} along with compressibility values
obtained from both other SPC/E simulations and experiment.

\begin{figure}
\includegraphics[width=\linewidth]{new_isothermalN}
\caption{Compressibility of SPC/E water}
\label{fig:compWater}
\end{figure}

Isothermal compressibility values calculated using the number density
(Eq. \ref{eq:BMN}) expression are in good agreement with experimental
and previous simulation work throughout the 1 -- 1000 atm pressure
regime.  Compressibilities computed using the Hull volume, however,
deviate dramatically from the experimental values at low applied
pressures.  The reason for this deviation is quite simple: at low
applied pressures, the liquid is in equilibrium with a vapor phase,
and it is entirely possible for one (or a few) molecules to drift away
from the liquid cluster (see Fig. \ref{fig:coneOfShame}).  At low
pressures, the restoring forces on the facets are very gentle, and
this means that the hulls often take on relatively distorted
geometries which include large volumes of empty space.

\begin{figure}
\includegraphics[width=\linewidth]{coneOfShame}
\caption{At low pressures, the liquid is in equilibrium with the vapor
  phase, and isolated molecules can detach from the liquid droplet.
  This is expected behavior, but the volume of the convex hull
  includes large regions of empty space. For this reason,
  compressibilities are computed using local number densities rather
  than hull volumes.}
\label{fig:coneOfShame}
\end{figure}

At higher pressures, the equilibrium strongly favors the liquid phase,
and the hull geometries are much more compact.  Because of the
liquid-vapor effect on the convex hull, the regional number density
approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the
compressibility.

In both the traditional compressibility formula (Eq. \ref{eq:BM}) and
the number density version (Eq. \ref{eq:BMN}), multiple simulations at
different pressures must be done to compute the first derivatives.  It
is also possible to compute the compressibility using the fluctuation
dissipation theorem using either fluctuations in the
volume,\cite{Debenedetti1986}
\begin{equation}
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle
    V \right \rangle ^{2}}{V \, k_{B} \, T},
\label{eq:BMVfluct}
\end{equation}
or, equivalently, fluctuations in the number of molecules within the
fixed region,
\begin{equation}
\kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle
    N \right \rangle ^{2}}{N \, k_{B} \, T}.
\label{eq:BMNfluct}
\end{equation}
Thus, the compressibility of each simulation can be calculated
entirely independently from other trajectories.  Compressibility
calculations that rely on the hull volume will still suffer the
effects of the empty space due to the vapor phase; for this reason, we
recommend using the number density (Eq. \ref{eq:BMN}) or number
density fluctuations (Eq. \ref{eq:BMNfluct}) for computing
compressibilities.  We obtained the results in
Fig. \ref{fig:compWater} using a sampling radius that was
approximately 80\% of the mean distance between the center of mass of
the cluster and the hull atoms.  This ratio of sampling radius to
average hull radius excludes the problematic vapor phase on the
outside of the cluster while including enough of the liquid phase to
avoid poor statistics due to fluctuating local densities.

A comparison of the oxygen-oxygen radial distribution functions for
SPC/E water simulated using both the Langevin Hull and more
traditional periodic boundary methods -- both at 1 atm and 300K --
reveals an understructuring of water in the Langevin Hull that
manifests as a slight broadening of the solvation shells.  This effect
may be due to the introduction of a liquid-vapor interface in the
Langevin Hull simulations (an interface which is missing in most
periodic simulations of bulk water).  Vapor-phase molecules contribute
a small but nearly flat portion of the radial distribution function.

\subsection{Molecular orientation distribution at cluster boundary}

In order for a non-periodic boundary method to be widely applicable,
it must be constructed in such a way that they allow a finite system
to replicate the properties of the bulk. Early non-periodic simulation
methods (e.g. hydrophobic boundary potentials) induced spurious
orientational correlations deep within the simulated
system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
fixing and characterizing the effects of artificial boundaries
including methods which fix the orientations of a set of edge
molecules.\cite{Warshel1978,King1989}

As described above, the Langevin Hull does not require that the
orientation of molecules be fixed, nor does it utilize an explicitly
hydrophobic boundary, or orientational or radial constraints.
Therefore, the orientational correlations of the molecules in water
clusters are of particular interest in testing this method.  Ideally,
the water molecules on the surfaces of the clusters will have enough
mobility into and out of the center of the cluster to maintain
bulk-like orientational distribution in the absence of orientational
and radial constraints.  However, since the number of hydrogen bonding
partners available to molecules on the exterior are limited, it is
likely that there will be an effective hydrophobicity of the hull.

To determine the extent of these effects, 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.  The
orientational angle of a water molecule is described by
\begin{equation} 
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
\end{equation}
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}.  Bulk-like
distributions will result in $\langle \cos \theta \rangle$ values
close to zero.  If the hull exhibits an overabundance of
externally-oriented oxygen sites, the average orientation will be
negative, while dangling hydrogen sites will result in positive
average orientations.

Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
for molecules in the interior of the cluster (squares) and for
molecules included in the convex hull (circles).
\begin{figure}
\includegraphics[width=\linewidth]{pAngle}
\caption{Distribution of $\cos{\theta}$ values for molecules on the
  interior of the cluster (squares) and for those participating in the
  convex hull (circles) at a variety of pressures.  The Langevin Hull
  exhibits minor dewetting behavior with exposed oxygen sites on the
  hull water molecules.  The orientational preference for exposed
  oxygen appears to be independent of applied pressure. }
\label{fig:pAngle}
\end{figure}

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 oxygen directed toward the exterior of the cluster,
forming dangling hydrogen bond acceptor sites.

The orientational preference exhibited by water molecules on the hull
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 near the cluster surface.

Previous molecular dynamics simulations of SPC/E liquid / vapor
interfaces using periodic boundary conditions have shown that
molecules on the liquid side of interface favor a similar orientation
where oxygen is directed away from the bulk.\cite{Taylor1996} These
simulations had well-defined liquid and vapor phase regions
equilibrium and it was observed that {\it vapor} molecules generally
had one hydrogen protruding from the surface, forming a dangling
hydrogen bond donor. Our water clusters do not have a true vapor
region, but rather a few transient molecules that leave the liquid
droplet (and which return to the droplet relatively quickly).
Although we cannot obtain an orientational preference of vapor phase
molecules in a Langevin Hull simulation, but we do agree with previous
estimates of the orientation of {\it liquid phase} molecules at the
interface.

\subsection{Heterogeneous nanoparticle / water mixtures}

To further test the method, we simulated gold nanoparticles ($r = 18$
\AA~, 1433 atoms) solvated by explicit SPC/E water clusters (5000
molecules) using a model for the gold / water interactions that has
been used by Dou {\it et. al.} for investigating the separation of
water films near hot metal surfaces.\cite{ISI:000167766600035} The
Langevin Hull was used to sample pressures of 1, 2, 5, 10, 20, 50, 100
and 200 atm, while all simulations were done at a temperature of 300
K.  At these temperatures and pressures, there is no observed
separation of the water film from the surface.

In Fig. \ref{fig:RhoR} we show the density of water and gold as a
function of the distance from the center of the nanoparticle.  Higher
applied pressures appear to destroy structural correlations in the
outermost monolayer of the gold nanoparticle as well as in the water
at the near the metal / water interface.  Simulations at increased
pressures exhibit significant overlap of the gold and water densities,
indicating a less well-defined interfacial surface.

\begin{figure}
\includegraphics[width=\linewidth]{RhoR}
\caption{Density profiles of gold and water at the nanoparticle
  surface. Each curve has been normalized by the average density in
  the bulk-like region available to the corresponding material.
  Higher applied pressures de-structure both the gold nanoparticle
  surface and water at the metal/water interface.}
\label{fig:RhoR}
\end{figure}

At even higher pressures (500 atm and above), problems with the metal
- water interaction potential became quite clear.  The model we are
using appears to have been parameterized for relatively low pressures;
it utilizes both shifted Morse and repulsive Morse potentials to model
the Au/O and Au/H interactions, respectively.  The repulsive wall of
the Morse potential does not diverge quickly enough at short distances
to prevent water from diffusing into the center of the gold
nanoparticles.  This behavior is likely not a realistic description of
the real physics of the situation.  A better model of the gold-water
adsorption behavior would require harder repulsive walls to prevent
this behavior.

\section{Discussion}
\label{sec:discussion}

The Langevin Hull samples the isobaric-isothermal ensemble for
non-periodic systems by coupling the system to a bath characterized by
pressure, temperature, and solvent viscosity.  This enables the
simulation of heterogeneous systems composed of materials with
significantly different compressibilities.  Because the boundary is
dynamically determined during the simulation and the molecules
interacting with the boundary can change, the method inflicts minimal
perturbations on the behavior of molecules at the edges of the
simulation.  Further work on this method will involve implicit
electrostatics at the boundary (which is missing in the current
implementation) as well as more sophisticated treatments of the
surface geometry (alpha
shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight
Cocone\cite{Dey:2003ts}). The non-convex hull geometries are
significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull
($\mathcal{O}(N \log N)$), but would enable the use of hull volumes
directly in computing the compressibility of the sample.

\section*{Appendix A: Computing Convex Hulls on Parallel Computers} 

In order to use the Langevin Hull for simulations on parallel
computers, one of the more difficult tasks is to compute the bounding
surface, facets, and resistance tensors when the individual processors
have incomplete information about the entire system's topology.  Most
parallel decomposition methods assign primary responsibility for the
motion of an atomic site to a single processor, and we can exploit
this to efficiently compute the convex hull for the entire system.

The basic idea involves splitting the point cloud into
spatially-overlapping subsets and computing the convex hulls for each
of the subsets.  The points on the convex hull of the entire system
are all present on at least one of the subset hulls. The algorithm
works as follows:
\begin{enumerate}
\item Each processor computes the convex hull for its own atomic sites
  (left panel in Fig. \ref{fig:parallel}).
\item The Hull vertices from each processor are communicated to all of
  the processors, and each processor assembles a complete list of hull
  sites (this is much smaller than the original number of points in
  the point cloud).
\item Each processor computes the global convex hull (right panel in
  Fig. \ref{fig:parallel}) using only those points that are the union
  of sites gathered from all of the subset hulls.  Delaunay
  triangulation is then done to obtain the facets of the global hull.
\end{enumerate}

\begin{figure}
\includegraphics[width=\linewidth]{parallel}
\caption{When the sites are distributed among many nodes for parallel
  computation, the processors first compute the convex hulls for their
  own sites (dashed lines in left panel). The positions of the sites
  that make up the subset hulls are then communicated to all
  processors (middle panel).  The convex hull of the system (solid line in 
  right panel) is the convex hull of the points on the union of the subset 
  hulls.}
\label{fig:parallel}
\end{figure}

The individual hull operations scale with
$\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total
number of sites, and $p$ is the number of processors.  These local
hull operations create a set of $p$ hulls, each with approximately
$\frac{n}{3pr}$ sites for a cluster of radius $r$. The worst-case
communication cost for using a ``gather'' operation to distribute this
information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n
  \beta (p-1)}{3 r p^2})$, while the final computation of the system
hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$.

For a large number of atoms on a moderately parallel machine, the
total costs are dominated by the computations of the individual hulls,
and communication of these hulls to create the Langevin Hull sees roughly
linear speed-up with increasing processor counts.

\section*{Acknowledgments}
Support for this project was provided by the
National Science Foundation under grant CHE-0848243. Computational
time was provided by the Center for Research Computing (CRC) at the
University of Notre Dame.  

Molecular graphics images were produced using the UCSF Chimera package from 
the Resource for Biocomputing, Visualization, and Informatics at the 
University of California, San Francisco (supported by NIH P41 RR001081).
\newpage

\bibliography{langevinHull}

\end{doublespace}
\end{document}
Revision:	3723
Committed:	Fri Feb 18 14:08:15 2011 UTC (13 years, 4 months ago) by gezelter
Content type:	application/x-tex
File size:	43474 byte(s)
Log Message:	final version