ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/langevinHull/langevinHull.tex
Revision: 3716
Committed: Thu Jan 27 00:00:20 2011 UTC (14 years, 7 months ago) by kstocke1
Content type: application/x-tex
File size: 43472 byte(s)
Log Message:
Made changes from reviews.

File Contents

# Content
1 \documentclass[11pt]{article}
2 \usepackage{amsmath}
3 \usepackage{amssymb}
4 \usepackage{setspace}
5 \usepackage{endfloat}
6 \usepackage{caption}
7 \usepackage{graphicx}
8 \usepackage{multirow}
9 \usepackage[square, comma, sort&compress]{natbib}
10 \usepackage{url}
11 \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
12 \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
13 9.0in \textwidth 6.5in \brokenpenalty=10000
14
15 % double space list of tables and figures
16 %\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
17 \setlength{\abovecaptionskip}{20 pt}
18 \setlength{\belowcaptionskip}{30 pt}
19
20 \bibpunct{}{}{,}{s}{}{;}
21 \bibliographystyle{achemso}
22
23 \begin{document}
24
25 \title{The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems}
26
27 \author{Charles F. Vardeman II, Kelsey M. Stocker, and J. Daniel
28 Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
29 Department of Chemistry and Biochemistry,\\
30 University of Notre Dame\\
31 Notre Dame, Indiana 46556}
32
33 \date{\today}
34
35 \maketitle
36
37 \begin{doublespace}
38
39 \begin{abstract}
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 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
58
59 %\narrowtext
60
61 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62 % BODY OF TEXT
63 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64
65
66 \section{Introduction}
67
68 The most common molecular dynamics methods for sampling configurations
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 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
92 are incompressible relative to other portions. Examples include
93 simulations of metallic nanoparticles in liquid environments, proteins
94 at ice / water interfaces, as well as other heterogeneous or
95 interfacial environments. In these cases, the affine transform of
96 atomic coordinates will either cause numerical instability when the
97 sites in the incompressible medium collide with each other, or will
98 lead to inefficient sampling of system volumes if the barostat is set
99 slow enough to avoid the instabilities in the incompressible region.
100
101 \begin{figure}
102 \includegraphics[width=\linewidth]{AffineScale2}
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
113 One may also wish to avoid affine transform periodic boundary methods
114 to simulate {\it explicitly non-periodic systems} under constant
115 pressure conditions. The use of periodic boxes to enforce a system
116 volume 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 boxes
119 solvating a protein under periodic boundary conditions are quite
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 {\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 structures.
130
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
134 nearly-constant {\it volume} conditions while maintaining bulk-like
135 behavior. Berkowitz and McCammon introduced a stochastic (Langevin)
136 boundary layer inside a region of fixed molecules which effectively
137 enforces constant temperature and volume (NVT)
138 conditions.\cite{Berkowitz1982} In this approach, the stochastic and
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.
143
144 The electrostatic and dispersive behavior near the boundary has long
145 been a cause for concern when performing simulations of explicitly
146 non-periodic systems. Early work led to the surface constrained soft
147 sphere dipole model (SCSSD)\cite{Warshel1978} in which the surface
148 molecules are fixed in a random orientation representative of the bulk
149 solvent structural properties. Belch {\it et al.}\cite{Belch1985}
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 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.
170
171 \subsection*{Restraining Potentials}
172 Restraining {\it potentials} introduce repulsive potentials at the
173 surface of a sphere or other geometry. The solute and any explicit
174 solvent are therefore restrained inside the range defined by the
175 external potential. Often the potentials include a weak short-range
176 attraction to maintain the correct density at the boundary. Beglov
177 and Roux have also introduced a restraining boundary potential which
178 relaxes dynamically depending on the solute geometry and the force the
179 explicit system exerts on the shell.\cite{Beglov:1995fk}
180
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: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
188 bulk solvent. Although this approach is appealing and has physical
189 motivation, nanoparticles do not deform far from their original
190 geometries even at temperatures which vaporize the nearby solvent. For
191 the systems like this, the flexible boundary model will be nearly
192 identical to a fixed-volume restraining potential.
193
194 \subsection*{Hull methods}
195 The approach of Kohanoff, Caro, and Finnis is the most promising of
196 the methods for introducing both constant pressure and temperature
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 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
205 bounding surface surrounding the outermost atoms in the simulated
206 system. This is not the only possible triangulated outer surface, but
207 guarantees that all of the random forces point inward towards the
208 cluster.
209
210 In the following sections, we extend and generalize the approach of
211 Kohanoff, Caro, and Finnis. The new method, which we are calling the
212 ``Langevin Hull'' applies the external pressure, Langevin drag, and
213 random forces on the {\it facets of the hull} instead of the atomic
214 sites comprising the vertices of the hull. This allows us to decouple
215 the external pressure contribution from the drag and random force.
216 The methodology is introduced in section \ref{sec:meth}, tests on
217 crystalline nanoparticles, liquid clusters, and heterogeneous mixtures
218 are detailed in section \ref{sec:tests}. Section \ref{sec:discussion}
219 summarizes our findings.
220
221 \section{Methodology}
222 \label{sec:meth}
223
224 The Langevin Hull uses an external bath at a fixed constant pressure
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
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
240 eliminates all facets on the interior of the point cloud, leaving only
241 those exposed to the bath. Sites on the convex hull are dynamic; as
242 molecules re-enter the cluster, all interactions between atoms on that
243 molecule and the external bath are removed. Since the edge is
244 determined dynamically as the simulation progresses, no {\it a priori}
245 geometry is defined. The pressure and temperature bath interacts only
246 with the atoms on the edge and not with atoms interior to the
247 simulation.
248
249 \begin{figure}
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.}
255 \label{fig:hullSample}
256 \end{figure}
257
258 Atomic sites in the interior of the simulation move under standard
259 Newtonian dynamics,
260 \begin{equation}
261 m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U,
262 \label{eq:Newton}
263 \end{equation}
264 where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the
265 instantaneous velocity of site $i$ at time $t$, and $U$ is the total
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}$:
270 \begin{equation}
271 m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}.
272 \end{equation}
273
274 The external bath interacts indirectly with the atomic sites through
275 the intermediary of the hull facets. Since each vertex (or atom)
276 provides one corner of a triangular facet, the force on the facets are
277 divided equally to each vertex. However, each vertex can participate
278 in multiple facets, so the resultant force is a sum over all facets
279 $f$ containing vertex $i$:
280 \begin{equation}
281 {\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\
282 } f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf
283 F}_f^{\mathrm ext}
284 \end{equation}
285
286 The external pressure bath applies a force to the facets of the convex
287 hull in direct proportion to the area of the facet, while the thermal
288 coupling depends on the solvent temperature, viscosity and the size
289 and shape of each facet. The thermal interactions are expressed as a
290 standard Langevin description of the forces,
291 \begin{equation}
292 \begin{array}{rclclcl}
293 {\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\
294 & = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t)
295 \end{array}
296 \end{equation}
297 Here, $A_f$ and $\hat{n}_f$ are the area and (outward-facing) normal
298 vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is the
299 velocity of the facet centroid,
300 \begin{equation}
301 {\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i,
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 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}
309 \left< {\mathbf R}_f(t) \right> & = & 0 \\
310 \left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\
311 \Xi_f(t)\delta(t-t^\prime).
312 \label{eq:randomForce}
313 \end{eqnarray}
314
315 Once the resistance tensor is known for a given facet, a stochastic
316 vector that has the properties in Eq. (\ref{eq:randomForce}) can be
317 calculated efficiently by carrying out a Cholesky decomposition to
318 obtain the square root matrix of the resistance tensor,
319 \begin{equation}
320 \Xi_f = {\bf S} {\bf S}^{T},
321 \label{eq:Cholesky}
322 \end{equation}
323 where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
324 vector with the statistics required for the random force can then be
325 obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which
326 has elements chosen from a Gaussian distribution, such that:
327 \begin{equation}
328 \langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
329 {\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
330 \end{equation}
331 where $\delta t$ is the timestep in use during the simulation. The
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_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
339 tensors for rigid bodies has been detailed
340 elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun:2008fk}
341 but the standard approach involving bead approximations would be
342 prohibitively expensive if it were recomputed at each step in a
343 molecular dynamics simulation.
344
345 Instead, we are utilizing an approximate resistance tensor obtained by
346 first constructing the Oseen tensor for the interaction of the
347 centroid of the facet ($f$) with each of the subfacets $\ell=1,2,3$,
348 \begin{equation}
349 T_{\ell f}=\frac{A_\ell}{8\pi\eta R_{\ell f}}\left(I +
350 \frac{\mathbf{R}_{\ell f}\mathbf{R}_{\ell f}^T}{R_{\ell f}^2}\right)
351 \end{equation}
352 Here, $A_\ell$ is the area of subfacet $\ell$ which is a triangle
353 containing two of the vertices of the facet along with the centroid.
354 $\mathbf{R}_{\ell f}$ is the vector between the centroid of facet $f$
355 and the centroid of sub-facet $\ell$, and $I$ is the ($3 \times 3$)
356 identity matrix. $\eta$ is the viscosity of the external bath.
357
358 \begin{figure}
359 \includegraphics[width=\linewidth]{hydro}
360 \caption{The resistance tensor $\Xi$ for a facet comprising sites $i$,
361 $j$, and $k$ is constructed using Oseen tensor contributions between
362 the centoid of the facet $f$ and each of the sub-facets ($i,f,j$),
363 ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets are
364 located at $1$, $2$, and $3$, and the area of each sub-facet is
365 easily computed using half the cross product of two of the edges.}
366 \label{hydro}
367 \end{figure}
368
369 The tensors for each of the sub-facets are added together, and the
370 resulting matrix is inverted to give a $3 \times 3$ resistance tensor
371 for translations of the triangular facet,
372 \begin{equation}
373 \Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}.
374 \end{equation}
375 Note that this treatment ignores rotations (and
376 translational-rotational coupling) of the facet. In compact systems,
377 the facets stay relatively fixed in orientation between
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,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.
385 \item Delaunay triangulation is carried out using the current atomic
386 configuration.
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 ($-\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.
394 \item[c.] Facet drag ($-\Xi_f(t) \mathbf{v}_f(t)$) forces are
395 computed.
396 \item[d.] Random forces ($\mathbf{R}_f(t)$) are computed using the
397 resistance tensor and the temperature ($T$) of the bath.
398 \end{itemize}
399 \item The facet forces are divided equally among the vertex atoms.
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{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
407 to parallelize on distributed memory machines (see Appendix A).
408
409 \section{Tests \& Applications}
410 \label{sec:tests}
411
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 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}.
419 \label{eq:BM}
420 \end{equation}
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
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
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}.
432 \label{eq:BMN}
433 \end{equation}
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 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
442 (Eq. \ref{eq:BM}), but this has implications at lower pressures (which
443 are explored in detail in the section on water droplets).
444
445 The metallic force field in use for the gold nanoparticles is the
446 quantum Sutton-Chen (QSC) model.\cite{PhysRevB.59.3527} In all
447 simulations involving point charges, we utilized damped shifted-force
448 (DSF) electrostatics\cite{Fennell06} which is a variant of the Wolf
449 summation\cite{wolf:8254} that has been shown to provide good forces
450 and torques on molecular models for water in a computationally
451 efficient manner.\cite{Fennell06} The damping parameter ($\alpha$) was
452 set to 0.18 \AA$^{-1}$, and the cutoff radius was set to 12 \AA. The
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{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 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
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] ,
474 \end{equation}
475 where $V^{pair}_{ij}$ and $\rho_{i}$ are given by
476 \begin{equation}
477 \label{eq:SCP2}
478 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}}.
479 \end{equation}
480 $V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for
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
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
488 tuned to various experimental properties such as the density, cohesive
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,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 crystal have yielded values
497 of 175.53 GPa.\cite{QSC2} Using the same force field, we have performed
498 a series of 1 ns simulations on gold nanoparticles of three different radii under 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} Polycrystalline gold has a reported bulk modulus of 220 GPa. 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 approach the polycrystalline bulk modulus value while larger nanoparticles approach the single crystal value. As nanoparticle size decreases, the bulk modulus becomes larger and the nanoparticle is less compressible. 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.
499
500 We measure 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}}$.
501
502 \begin{figure}
503 \includegraphics[width=\linewidth]{stacked}
504 \caption{The response of the internal pressure and temperature of gold
505 nanoparticles when first placed in the Langevin Hull
506 ($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting
507 from initial conditions that were far from the bath pressure and
508 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).}
509 \label{fig:pressureResponse}
510 \end{figure}
511
512 We note that the Langevin Hull produces rapidly-converging behavior
513 for structures that are started far from equilibrium. In
514 Fig. \ref{fig:pressureResponse} we show how the pressure and
515 temperature respond to the Langevin Hull for nanoparticles that were
516 initialized far from the target pressure and temperature. As
517 expected, the rate at which thermal equilibrium is achieved depends on
518 the total surface area of the cluster exposed to the bath as well as
519 the bath viscosity. Pressure that is applied suddenly to a cluster
520 can excite breathing vibrations, but these rapidly damp out (on time
521 scales of 30 -- 50 ps).
522
523 \subsection{Compressibility of SPC/E water clusters}
524
525 Prior molecular dynamics simulations on SPC/E water (both in
526 NVT~\cite{Glattli2002} and NPT~\cite{Motakabbir1990, Pi2009}
527 ensembles) have yielded values for the isothermal compressibility that
528 agree well with experiment.\cite{Fine1973} The results of two
529 different approaches for computing the isothermal compressibility from
530 Langevin Hull simulations for pressures between 1 and 3000 atm are
531 shown in Fig. \ref{fig:compWater} along with compressibility values
532 obtained from both other SPC/E simulations and experiment.
533
534 \begin{figure}
535 \includegraphics[width=\linewidth]{new_isothermalN}
536 \caption{Compressibility of SPC/E water}
537 \label{fig:compWater}
538 \end{figure}
539
540 Isothermal compressibility values calculated using the number density
541 (Eq. \ref{eq:BMN}) expression are in good agreement with experimental
542 and previous simulation work throughout the 1 -- 1000 atm pressure
543 regime. Compressibilities computed using the Hull volume, however,
544 deviate dramatically from the experimental values at low applied
545 pressures. The reason for this deviation is quite simple: at low
546 applied pressures, the liquid is in equilibrium with a vapor phase,
547 and it is entirely possible for one (or a few) molecules to drift away
548 from the liquid cluster (see Fig. \ref{fig:coneOfShame}). At low
549 pressures, the restoring forces on the facets are very gentle, and
550 this means that the hulls often take on relatively distorted
551 geometries which include large volumes of empty space.
552
553 \begin{figure}
554 \includegraphics[width=\linewidth]{coneOfShame}
555 \caption{At low pressures, the liquid is in equilibrium with the vapor
556 phase, and isolated molecules can detach from the liquid droplet.
557 This is expected behavior, but the volume of the convex hull
558 includes large regions of empty space. For this reason,
559 compressibilities are computed using local number densities rather
560 than hull volumes.}
561 \label{fig:coneOfShame}
562 \end{figure}
563
564 At higher pressures, the equilibrium strongly favors the liquid phase,
565 and the hull geometries are much more compact. Because of the
566 liquid-vapor effect on the convex hull, the regional number density
567 approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the
568 compressibility.
569
570 In both the traditional compressibility formula (Eq. \ref{eq:BM}) and
571 the number density version (Eq. \ref{eq:BMN}), multiple simulations at
572 different pressures must be done to compute the first derivatives. It
573 is also possible to compute the compressibility using the fluctuation
574 dissipation theorem using either fluctuations in the
575 volume,\cite{Debenedetti1986}
576 \begin{equation}
577 \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle
578 V \right \rangle ^{2}}{V \, k_{B} \, T},
579 \label{eq:BMVfluct}
580 \end{equation}
581 or, equivalently, fluctuations in the number of molecules within the
582 fixed region,
583 \begin{equation}
584 \kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle
585 N \right \rangle ^{2}}{N \, k_{B} \, T}.
586 \label{eq:BMNfluct}
587 \end{equation}
588 Thus, the compressibility of each simulation can be calculated
589 entirely independently from other trajectories. Compressibility
590 calculations that rely on the hull volume will still suffer the
591 effects of the empty space due to the vapor phase; for this reason, we
592 recommend using the number density (Eq. \ref{eq:BMN}) or number
593 density fluctuations (Eq. \ref{eq:BMNfluct}) for computing
594 compressibilities. We achieved the best results using a sampling radius approximately 80\% of the cluster radius. This ratio of sampling radius to cluster 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.
595
596 A comparison of the oxygen-oxygen radial distribution functions for SPC/E water simulated using the Langevin Hull and bulk SPC/E using periodic boundary conditions -- both at 1 atm and 300K -- reveals a slight understructuring of water in the Langevin Hull that manifests as a minor broadening of the solvation shells. This effect may be related to the introduction of surface tension around the entire cluster, an effect absent in bulk systems. As a result, molecules on the hull may experience an increased inward force, slightly compressing the solvation shell structure.
597
598 \subsection{Molecular orientation distribution at cluster boundary}
599
600 In order for a non-periodic boundary method to be widely applicable,
601 it must be constructed in such a way that they allow a finite system
602 to replicate the properties of the bulk. Early non-periodic simulation
603 methods (e.g. hydrophobic boundary potentials) induced spurious
604 orientational correlations deep within the simulated
605 system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
606 fixing and characterizing the effects of artificial boundaries
607 including methods which fix the orientations of a set of edge
608 molecules.\cite{Warshel1978,King1989}
609
610 As described above, the Langevin Hull does not require that the
611 orientation of molecules be fixed, nor does it utilize an explicitly
612 hydrophobic boundary, or orientational or radial constraints.
613 Therefore, the orientational correlations of the molecules in water
614 clusters are of particular interest in testing this method. Ideally,
615 the water molecules on the surfaces of the clusters will have enough
616 mobility into and out of the center of the cluster to maintain
617 bulk-like orientational distribution in the absence of orientational
618 and radial constraints. However, since the number of hydrogen bonding
619 partners available to molecules on the exterior are limited, it is
620 likely that there will be an effective hydrophobicity of the hull.
621
622 To determine the extent of these effects, we examined the
623 orientations exhibited by SPC/E water in a cluster of 1372
624 molecules at 300 K and at pressures ranging from 1 -- 1000 atm. The
625 orientational angle of a water molecule is described by
626 \begin{equation}
627 \cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
628 \end{equation}
629 where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of
630 mass and the cluster center of mass, and $\vec{\mu}_{i}$ is the vector
631 bisecting the H-O-H angle of molecule {\it i}. Bulk-like
632 distributions will result in $\langle \cos \theta \rangle$ values
633 close to zero. If the hull exhibits an overabundance of
634 externally-oriented oxygen sites, the average orientation will be
635 negative, while dangling hydrogen sites will result in positive
636 average orientations.
637
638 Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
639 for molecules in the interior of the cluster (squares) and for
640 molecules included in the convex hull (circles).
641 \begin{figure}
642 \includegraphics[width=\linewidth]{pAngle}
643 \caption{Distribution of $\cos{\theta}$ values for molecules on the
644 interior of the cluster (squares) and for those participating in the
645 convex hull (circles) at a variety of pressures. The Langevin Hull
646 exhibits minor dewetting behavior with exposed oxygen sites on the
647 hull water molecules. The orientational preference for exposed
648 oxygen appears to be independent of applied pressure. }
649 \label{fig:pAngle}
650 \end{figure}
651
652 As expected, interior molecules (those not included in the convex
653 hull) maintain a bulk-like structure with a uniform distribution of
654 orientations. Molecules included in the convex hull show a slight
655 preference for values of $\cos{\theta} < 0.$ These values correspond
656 to molecules with oxygen directed toward the exterior of the cluster,
657 forming dangling hydrogen bond acceptor sites.
658
659 The orientational preference exhibited by water molecules on the hull
660 is significantly weaker than the preference caused by an explicit
661 hydrophobic bounding potential. Additionally, the Langevin Hull does
662 not require that the orientation of any molecules be fixed in order to
663 maintain bulk-like structure, even near the cluster surface.
664
665 Previous molecular dynamics simulations of SPC/E liquid / vapor
666 interfaces using periodic boundary conditions have shown that
667 molecules on the liquid side of interface favor a similar orientation
668 where oxygen is directed away from the bulk.\cite{Taylor1996} These
669 simulations had well-defined liquid and vapor phase regions
670 equilibrium and it was observed that {\it vapor} molecules generally
671 had one hydrogen protruding from the surface, forming a dangling
672 hydrogen bond donor. Our water clusters do not have a true vapor
673 region, but rather a few transient molecules that leave the liquid
674 droplet (and which return to the droplet relatively quickly).
675 Although we cannot obtain an orientational preference of vapor phase
676 molecules in a Langevin Hull simulation, but we do agree with previous
677 estimates of the orientation of {\it liquid phase} molecules at the
678 interface.
679
680 \subsection{Heterogeneous nanoparticle / water mixtures}
681
682 To further test the method, we simulated gold nanoparticles ($r = 18$
683 \AA) solvated by explicit SPC/E water clusters using a model for the
684 gold / water interactions that has been used by Dou {\it et. al.} for
685 investigating the separation of water films near hot metal
686 surfaces.\cite{ISI:000167766600035} The Langevin Hull was used to
687 sample pressures of 1, 2, 5, 10, 20, 50, 100 and 200 atm, while all
688 simulations were done at a temperature of 300 K. At these
689 temperatures and pressures, there is no observed separation of the
690 water film from the surface.
691
692 In Fig. \ref{fig:RhoR} we show the density of water and gold as a
693 function of the distance from the center of the nanoparticle. Higher
694 applied pressures appear to destroy structural correlations in the
695 outermost monolayer of the gold nanoparticle as well as in the water
696 at the near the metal / water interface. Simulations at increased
697 pressures exhibit significant overlap of the gold and water densities,
698 indicating a less well-defined interfacial surface.
699
700 \begin{figure}
701 \includegraphics[width=\linewidth]{RhoR}
702 \caption{Density profiles of gold and water at the nanoparticle
703 surface. Each curve has been normalized by the average density in
704 the bulk-like region available to the corresponding material. Higher applied pressures
705 de-structure both the gold nanoparticle surface and water at the
706 metal/water interface.}
707 \label{fig:RhoR}
708 \end{figure}
709
710 At even higher pressures (500 atm and above), problems with the metal
711 - water interaction potential became quite clear. The model we are
712 using appears to have been parameterized for relatively low pressures;
713 it utilizes both shifted Morse and repulsive Morse potentials to model
714 the Au/O and Au/H interactions, respectively. The repulsive wall of
715 the Morse potential does not diverge quickly enough at short distances
716 to prevent water from diffusing into the center of the gold
717 nanoparticles. This behavior is likely not a realistic description of
718 the real physics of the situation. A better model of the gold-water
719 adsorption behavior appears to require harder repulsive walls to
720 prevent this behavior.
721
722 \section{Discussion}
723 \label{sec:discussion}
724
725 The Langevin Hull samples the isobaric-isothermal ensemble for
726 non-periodic systems by coupling the system to a bath characterized by
727 pressure, temperature, and solvent viscosity. This enables the
728 simulation of heterogeneous systems composed of materials with
729 significantly different compressibilities. Because the boundary is
730 dynamically determined during the simulation and the molecules
731 interacting with the boundary can change, the method inflicts minimal
732 perturbations on the behavior of molecules at the edges of the
733 simulation. Further work on this method will involve implicit
734 electrostatics at the boundary (which is missing in the current
735 implementation) as well as more sophisticated treatments of the
736 surface geometry (alpha
737 shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight
738 Cocone\cite{Dey:2003ts}). The non-convex hull geometries are
739 significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull
740 ($\mathcal{O}(N \log N)$), but would enable the use of hull volumes
741 directly in computing the compressibility of the sample.
742
743 \section*{Appendix A: Computing Convex Hulls on Parallel Computers}
744
745 In order to use the Langevin Hull for simulations on parallel
746 computers, one of the more difficult tasks is to compute the bounding
747 surface, facets, and resistance tensors when the individual processors
748 have incomplete information about the entire system's topology. Most
749 parallel decomposition methods assign primary responsibility for the
750 motion of an atomic site to a single processor, and we can exploit
751 this to efficiently compute the convex hull for the entire system.
752
753 The basic idea involves splitting the point cloud into
754 spatially-overlapping subsets and computing the convex hulls for each
755 of the subsets. The points on the convex hull of the entire system
756 are all present on at least one of the subset hulls. The algorithm
757 works as follows:
758 \begin{enumerate}
759 \item Each processor computes the convex hull for its own atomic sites
760 (left panel in Fig. \ref{fig:parallel}).
761 \item The Hull vertices from each processor are communicated to all of
762 the processors, and each processor assembles a complete list of hull
763 sites (this is much smaller than the original number of points in
764 the point cloud).
765 \item Each processor computes the global convex hull (right panel in
766 Fig. \ref{fig:parallel}) using only those points that are the union
767 of sites gathered from all of the subset hulls. Delaunay
768 triangulation is then done to obtain the facets of the global hull.
769 \end{enumerate}
770
771 \begin{figure}
772 \includegraphics[width=\linewidth]{parallel}
773 \caption{When the sites are distributed among many nodes for parallel
774 computation, the processors first compute the convex hulls for their
775 own sites (dashed lines in left panel). The positions of the sites
776 that make up the subset hulls are then communicated to all
777 processors (middle panel). The convex hull of the system (solid line in
778 right panel) is the convex hull of the points on the union of the subset
779 hulls.}
780 \label{fig:parallel}
781 \end{figure}
782
783 The individual hull operations scale with
784 $\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total
785 number of sites, and $p$ is the number of processors. These local
786 hull operations create a set of $p$ hulls, each with approximately
787 $\frac{n}{3pr}$ sites for a cluster of radius $r$. The worst-case
788 communication cost for using a ``gather'' operation to distribute this
789 information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n
790 \beta (p-1)}{3 r p^2})$, while the final computation of the system
791 hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$.
792
793 For a large number of atoms on a moderately parallel machine, the
794 total costs are dominated by the computations of the individual hulls,
795 and communication of these hulls to create the Langevin Hull sees roughly
796 linear speed-up with increasing processor counts.
797
798 \section*{Acknowledgments}
799 Support for this project was provided by the
800 National Science Foundation under grant CHE-0848243. Computational
801 time was provided by the Center for Research Computing (CRC) at the
802 University of Notre Dame.
803
804 Molecular graphics images were produced using the UCSF Chimera package from
805 the Resource for Biocomputing, Visualization, and Informatics at the
806 University of California, San Francisco (supported by NIH P41 RR001081).
807 \newpage
808
809 \bibliography{langevinHull}
810
811 \end{doublespace}
812 \end{document}