ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/langevinHull/langevinHull.tex
(Generate patch)

Comparing trunk/langevinHull/langevinHull.tex (file contents):
Revision 3663 by gezelter, Thu Oct 21 16:24:13 2010 UTC vs.
Revision 3723 by gezelter, Fri Feb 18 14:08:15 2011 UTC

# Line 17 | Line 17
17   \setlength{\abovecaptionskip}{20 pt}
18   \setlength{\belowcaptionskip}{30 pt}
19  
20 < \bibpunct{[}{]}{,}{s}{}{;}
21 < \bibliographystyle{aip}
20 > \bibpunct{}{}{,}{s}{}{;}
21 > \bibliographystyle{achemso}
22  
23   \begin{document}
24  
# Line 39 | Line 39 | Notre Dame, Indiana 46556}
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 objects in the system. Additionally, a Langevin
43 <  thermostat is applied to facets of the hull to mimic contact with an
44 <  external heat bath. This new method, the ``Langevin Hull'', performs
45 <  better than traditional affine transform methods for systems
46 <  containing heterogeneous mixtures of materials with different
47 <  compressibilities. It does not suffer from the edge effects of
48 <  boundary potential methods, and allows realistic treatment of both
49 <  external pressure and thermal conductivity to an implicit solvent.
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 systems are
46 >  problematic for traditional affine transform methods.  The Langevin
47 >  Hull does not suffer from the edge effects of boundary potential
48 >  methods, and allows realistic treatment of both external pressure
49 >  and thermal conductivity due to the presence of an implicit solvent.
50    We apply this method to several different systems including bare
51 <  nanoparticles, nanoparticles in an explicit solvent, as well as
52 <  clusters of liquid water and ice. The predicted mechanical and
53 <  thermal properties of these systems are in good agreement with
54 <  experimental data.
51 >  metal nanoparticles, nanoparticles in an explicit solvent, as well
52 >  as clusters of liquid water. The predicted mechanical properties of
53 >  these systems are in good agreement with experimental data and
54 >  previous simulation work.
55   \end{abstract}
56  
57   \newpage
# Line 66 | Line 66 | of an isobaric-isothermal (NPT) ensemble attempt to ma
66   \section{Introduction}
67  
68   The most common molecular dynamics methods for sampling configurations
69 < of an isobaric-isothermal (NPT) ensemble attempt to maintain a target
70 < pressure in a simulation by coupling the volume of the system to an
71 < extra degree of freedom, the {\it barostat}.  These methods require
72 < periodic boundary conditions, because when the instantaneous pressure
73 < in the system differs from the target pressure, the volume is
74 < typically reduced or expanded using {\it affine transforms} of the
75 < system geometry. An affine transform scales both the box lengths as
76 < well as the scaled particle positions (but not the sizes of the
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 coordinate
84 < transformation to adjust the box volume.  As long as the material in
85 < the simulation box is essentially a bulk-like liquid which has a
86 < relatively uniform compressibility, the standard affine transform
87 < approach provides an excellent way of adjusting the volume of the
88 < system and applying pressure directly via the interactions between
88 < atomic sites.
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 < The problem with this approach becomes apparent when the material
91 < being simulated is an inhomogeneous mixture in which portions of the
92 < simulation box are incompressible relative to other portions.
93 < Examples include simulations of metallic nanoparticles in liquid
94 < environments, proteins at interfaces, as well as other multi-phase or
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 lead to
98 < inefficient sampling of system volumes if the barostat is set slow
99 < enough to avoid the instabilities in the incompressible region.
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 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  
113   One may also wish to avoid affine transform periodic boundary methods
114   to simulate {\it explicitly non-periodic systems} under constant
115   pressure conditions. The use of periodic boxes to enforce a system
116 < volume either requires effective solute concentrations that are much
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 < There have been a number of other approaches to explicit
125 < non-periodicity that focus on constant or nearly-constant {\it volume}
126 < conditions while maintaining bulk-like behavior.  Berkowitz and
127 < McCammon introduced a stochastic (Langevin) boundary layer inside a
128 < region of fixed molecules which effectively enforces constant
129 < temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this
129 < approach, the stochastic and fixed regions were defined relative to a
130 < central atom.  Brooks and Karplus extended this method to include
131 < deformable stochastic boundaries.\cite{iii:6312} The stochastic
132 < boundary approach has been used widely for protein
133 < simulations. [CITATIONS NEEDED]
124 > {\it Total} protein concentrations in the cell are typically on the
125 > order of 160-310 mg/ml,\cite{Brown1991195} and individual proteins
126 > have concentrations orders of magnitude lower than this in the
127 > cellular environment. The effective concentrations of single proteins
128 > in simulations may have significant effects on the structure and
129 > dynamics of simulated systems.
130  
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.  King and Warshel introduced a surface
146 < constrained all-atom solvent (SCAAS) which included polarization
147 < effects of a fixed spherical boundary to mimic bulk-like behavior
148 < without periodic boundaries.\cite{king:3647} In the SCAAS model, a
149 < layer of fixed solvent molecules surrounds the solute and any explicit
150 < solvent, and this in turn is surrounded by a continuum dielectric.
151 < MORE HERE.  WHAT DID THEY FIND?
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 developed a boundary model in which the hard sphere
166 < boundary has a radius that varies with the instantaneous configuration
167 < of the solute (and solvent) molecules.\cite{beglov:9050} This model
168 < contains a clear pressure and surface tension contribution to the free
169 < energy which XXX.
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 this potential.  Often the
175 < potentials include a weak short-range attraction to maintain the
176 < correct density at the boundary.  Beglov and Roux have also introduced
177 < a restraining boundary potential which relaxes dynamically depending
178 < on the solute geometry and the force the explicit system exerts on the
179 < shell.\cite{Beglov:1995fk}
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 flexible boundary model
182 < that uses a Lennard-Jones potential between the solvent molecules and
183 < a boundary which is determined dynamically from the position of the
184 < nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This approach allows
185 < the confining potential to prevent solvent molecules from migrating
186 < too far from the solute surface, while providing a weak attractive
187 < force pulling the solvent molecules towards a fictitious bulk solvent.
188 < Although this approach is appealing and has physical motivation,
189 < nanoparticles do not deform far from their original geometries even at
190 < temperatures which vaporize the nearby solvent. For the systems like
191 < the one described, the flexible boundary model will be nearly
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 < force in a direction that is inward-facing relative to the facets of a
201 < closed bounding surface.  The statistical distribution of the random
200 > forces in a direction that is inward-facing relative to the facets of
201 > a closed bounding surface.  The statistical distribution of the random
202   forces are uniquely tied to the pressure in the external reservoir, so
203   the method can be shown to sample the isobaric-isothermal ensemble.
204   Kohanoff {\it et al.} used a Delaunay tessellation to generate a
# Line 187 | Line 210 | random forces on the facets of the {\it hull itself} i
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 facets of the {\it hull itself} instead of the
214 < atomic sites comprising the vertices of the hull.  This allows us to
215 < decouple the external pressure contribution from the drag and random
216 < force.  Section \ref{sec:meth}
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 < We have developed a new method which uses an external bath at a fixed
225 < constant pressure ($P$) and temperature ($T$).  This bath interacts
226 < only with the objects on the exterior hull of the system.  Defining
227 < the hull of the simulation is done in a manner similar to the approach
228 < of Kohanoff, Caro and Finnis.\cite{Kohanoff:2005qm} That is, any
229 < instantaneous configuration of the atoms in the system is considered
230 < as a point cloud in three dimensional space.  Delaunay triangulation
231 < is used to find all facets between coplanar neighbors.\cite{DELAUNAY}
232 < In highly symmetric point clouds, facets can contain many atoms, but
233 < in all but the most symmetric of cases the facets are simple triangles
234 < in 3-space that contain exactly three atoms.
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.  This eliminates all facets on the
240 < interior of the point cloud, leaving only those exposed to the
241 < bath. Sites on the convex hull are dynamic. As molecules re-enter the
242 < cluster, all interactions between atoms on that molecule and the
243 < external bath are removed.  Since the edge is determined dynamically
244 < as the simulation progresses, no {\it a priori} geometry is
245 < defined. The pressure and temperature bath interacts {\it directly}
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  
221 –
249   \begin{figure}
250 < \includegraphics[width=\linewidth]{hullSample}
250 > \includegraphics[width=\linewidth]{solvatedNano}
251   \caption{The external temperature and pressure bath interacts only
252    with those atoms on the convex hull (grey surface).  The hull is
253 <  computed dynamically at each time step, and molecules dynamically
254 <  move between the interior (Newtonian) region and the Langevin hull.}
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 <
232 < Atomic sites in the interior of the point cloud move under standard
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,
# Line 240 | Line 266 | equation of motion is modified with an external force,
266   potential energy.  For atoms on the exterior of the cluster
267   (i.e. those that occupy one of the vertices of the convex hull), the
268   equation of motion is modified with an external force, ${\mathbf
269 <  F}_i^{\mathrm ext}$,
269 >  F}_i^{\mathrm ext}$:
270   \begin{equation}
271   m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}.
272   \end{equation}
273  
274 < The external bath interacts directly with the facets of the convex
275 < hull. Since each vertex (or atom) provides one corner of a triangular
276 < facet, the force on the facets are divided equally to each vertex.
277 < However, each vertex can participate in multiple facets, so the resultant
278 < force is a sum over all facets $f$ containing vertex $i$:
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
# Line 267 | Line 294 | Here, $A_f$ and $\hat{n}_f$ are the area and normal ve
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 normal vectors for facet
298 < $f$, respectively.  ${\mathbf v}_f(t)$ is the velocity of the facet
299 < centroid,
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 fluid.  The resistance tensor is related to the
305 > viscosity of the bath.  The resistance tensor is related to the
306   fluctuations of the random force, $\mathbf{R}(t)$, by the
307   fluctuation-dissipation theorem,
308   \begin{eqnarray}
# Line 285 | Line 312 | Once the resistance tensor is known for a given facet
312   \label{eq:randomForce}
313   \end{eqnarray}
314  
315 < Once the resistance tensor is known for a given facet a stochastic
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 < done efficiently by carrying out a Cholesky decomposition to obtain
318 < the square root matrix of the resistance tensor,
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}
# Line 305 | Line 332 | Our treatment of the resistance tensor is approximate.
332   random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to
333   have the correct properties required by Eq. (\ref{eq:randomForce}).
334  
335 < Our treatment of the resistance tensor is approximate.  $\Xi$ for a
335 > Our treatment of the resistance tensor is approximate.  $\Xi_f$ for a
336   rigid triangular plate would normally be treated as a $6 \times 6$
337   tensor that includes translational and rotational drag as well as
338   translational-rotational coupling. The computation of resistance
# Line 315 | Line 342 | We are utilizing an approximate resistance tensor obta
342   prohibitively expensive if it were recomputed at each step in a
343   molecular dynamics simulation.
344  
345 < We are utilizing an approximate resistance tensor obtained by first
346 < constructing the Oseen tensor for the interaction of the centroid of
347 < the facet ($f$) with each of the subfacets $j$,
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_{jf}=\frac{A_j}{8\pi\eta R_{jf}}\left(I +
350 <  \frac{\mathbf{R}_{jf}\mathbf{R}_{jf}^T}{R_{jf}^2}\right)
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_j$ is the area of subfacet $j$ which is a triangle containing
353 < two of the vertices of the facet along with the centroid.
354 < $\mathbf{R}_{jf}$ is the vector between the centroid of facet $f$ and
355 < the centroid of sub-facet $j$, and $I$ is the ($3 \times 3$) identity
356 < matrix.  $\eta$ is the viscosity of the external bath.
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}
# Line 339 | Line 366 | The Oseen tensors for each of the sub-facets are added
366   \label{hydro}
367   \end{figure}
368  
369 < The Oseen tensors for each of the sub-facets are added together, and
370 < the resulting matrix is inverted to give a $3 \times 3$ resistance
371 < tensor for translations of the triangular facet,
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 explicitly ignores rotations (and
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,openmd} The Delaunay
382 < triangulation and computation of the convex hull are done using calls
383 < to the qhull library.\cite{Qhull} There is a moderate penalty for
384 < computing the convex hull at each step in the molecular dynamics
385 < simulation (HOW MUCH?), but the convex hull is remarkably easy to
386 < parallelize on distributed memory machines (see Appendix A).
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{SPCE} and 3) a heterogeneous
415 < mixture (gold nanoparticles in a water droplet). In each case, we have
368 < computed properties that depend on the external applied pressure.  Of
369 < particular interest for the single-phase systems is the bulk modulus,
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}.
# Line 375 | Line 421 | is not well-defined.  In order to compute the compress
421  
422   One problem with eliminating periodic boundary conditions and
423   simulation boxes is that the volume of a three-dimensional point cloud
424 < is not well-defined.  In order to compute the compressibility of a
424 > is not well-defined. In order to compute the compressibility of a
425   bulk material, we make an assumption that the number density, $\rho =
426 < \frac{N}{V}$, is uniform within some region of the cloud.  The
426 > \frac{N}{V}$, is uniform within some region of the point cloud. The
427   compressibility can then be expressed in terms of the average number
428   of particles in that region,
429   \begin{equation}
430 < \kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right
431 < )_{T}
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 10 \AA\ radius centered in
435 < the middle of the cluster. $N$ is the average number of molecules
434 > The region we used is a spherical volume of 20 \AA\ radius centered in
435 > the middle of the cluster with a roughly 25 \AA\ radius. $N$ is the average number of molecules
436   found within this region throughout a given simulation. The geometry
437 < and size of the region is arbitrary, and any bulk-like portion of the
438 < cluster can be used to compute the bulk modulus.
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 be taken as
441 < the system volume in the compressibility expression (Eq. \ref{eq:BM}),
442 < but this has implications at lower pressures (which are explored in
443 < detail in the section on water droplets).
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
# Line 407 | Line 453 | atoms and the SPC/E water molecules.\cite{ISI:00016776
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}
456 > \subsection{Bulk Modulus of gold nanoparticles}
457  
458 < The bulk modulus is well-known for gold, and it provides a good first
459 < test of how the method compares to other similar methods.  
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 <
467 < \begin{figure}
468 < \includegraphics[width=\linewidth]{pressure_tb}
469 < \caption{Pressure response is rapid (18 \AA gold nanoparticle), target
470 < pressure = 4 GPa}
420 < \label{pressureResponse}
421 < \end{figure}
422 <
423 < \begin{figure}
424 < \includegraphics[width=\linewidth]{temperature_tb}
425 < \caption{Temperature equilibration depends on surface area and bath
426 <  viscosity.  Target Temperature = 300K}
427 < \label{temperatureResponse}
428 < \end{figure}
429 <
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 < \kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial
473 <    P}\right)
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}$ and $D_{ii}$ set the appropriate overall energy
485 + scale, $c_i$ scales the attractive portion of the potential relative
486 + to the repulsive interaction and $\alpha_{ij}$ is a length parameter
487 + that assures a dimensionless form for $\rho$. These parameters are
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
498 + performed a series of 1 ns simulations on gold nanoparticles of three
499 + different radii: 20 \AA~ (1985 atoms), 30 \AA~ (6699 atoms), and 40
500 + \AA~ (15707 atoms) utilizing the Langevin Hull at a variety of applied
501 + pressures ranging from 0 -- 10 GPa.  For the 40 \AA~ radius
502 + nanoparticle we obtain a value of 177.55 GPa for the bulk modulus of
503 + gold, in close agreement with both previous simulations and the
504 + experimental bulk modulus reported for gold single
505 + crystals.\cite{Collard1991} The smaller gold nanoparticles (30 and 20
506 + \AA~ radii) have calculated bulk moduli of 215.58 and 208.86 GPa,
507 + respectively, indicating that smaller nanoparticles are somewhat
508 + stiffer (less compressible) than the larger nanoparticles.  This
509 + stiffening of the small nanoparticles may be related to their high
510 + degree of surface curvature, resulting in a lower coordination number
511 + of surface atoms relative to the the surface atoms in the 40 \AA~
512 + radius particle.
513 +
514 + We obtain a gold lattice constant of 4.051 \AA~ using the Langevin
515 + Hull at 1 atm, close to the experimentally-determined value for bulk
516 + gold and the value for gold simulated using the QSC potential and
517 + periodic boundary conditions (4.079 \AA~ and 4.088\AA~,
518 + respectively).\cite{QSC2} The slightly smaller calculated lattice
519 + constant is most likely due to the presence of surface tension in the
520 + non-periodic Langevin Hull cluster, an effect absent from a bulk
521 + simulation. The specific heat of a 40 \AA~ gold nanoparticle under the
522 + Langevin Hull at 1 atm is 24.914 $\mathrm {\frac{J}{mol \, K}}$, which
523 + compares very well with the experimental value of 25.42 $\mathrm
524 + {\frac{J}{mol \, K}}$.
525 +
526   \begin{figure}
527 < \includegraphics[width=\linewidth]{compress_tb}
528 < \caption{Isothermal Compressibility (18 \AA gold nanoparticle)}
529 < \label{temperatureResponse}
527 > \includegraphics[width=\linewidth]{stacked}
528 > \caption{The response of the internal pressure and temperature of gold
529 >  nanoparticles when first placed in the Langevin Hull
530 >  ($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting
531 >  from initial conditions that were far from the bath pressure and
532 >  temperature.  The pressure response is rapid (after the breathing mode oscillations in the nanoparticle die out), and the rate of thermal equilibration depends on both exposed surface area (top panel) and the viscosity of the bath (middle panel).}
533 > \label{fig:pressureResponse}
534   \end{figure}
535  
536 + We note that the Langevin Hull produces rapidly-converging behavior
537 + for structures that are started far from equilibrium.  In
538 + Fig. \ref{fig:pressureResponse} we show how the pressure and
539 + temperature respond to the Langevin Hull for nanoparticles that were
540 + initialized far from the target pressure and temperature.  As
541 + expected, the rate at which thermal equilibrium is achieved depends on
542 + the total surface area of the cluster exposed to the bath as well as
543 + the bath viscosity.  Pressure that is applied suddenly to a cluster
544 + can excite breathing vibrations, but these rapidly damp out (on time
545 + scales of 30 -- 50 ps).
546 +
547   \subsection{Compressibility of SPC/E water clusters}
548  
549   Prior molecular dynamics simulations on SPC/E water (both in
# Line 445 | Line 551 | Langevin Hull simulations for pressures between 1 and
551   ensembles) have yielded values for the isothermal compressibility that
552   agree well with experiment.\cite{Fine1973} The results of two
553   different approaches for computing the isothermal compressibility from
554 < Langevin Hull simulations for pressures between 1 and 6500 atm are
554 > Langevin Hull simulations for pressures between 1 and 3000 atm are
555   shown in Fig. \ref{fig:compWater} along with compressibility values
556   obtained from both other SPC/E simulations and experiment.
451 – Compressibility values from all references are for applied pressures
452 – within the range 1 - 1000 atm.
557  
558   \begin{figure}
559   \includegraphics[width=\linewidth]{new_isothermalN}
# Line 459 | Line 563 | and previous simulation work throughout the 1 - 1000 a
563  
564   Isothermal compressibility values calculated using the number density
565   (Eq. \ref{eq:BMN}) expression are in good agreement with experimental
566 < and previous simulation work throughout the 1 - 1000 atm pressure
566 > and previous simulation work throughout the 1 -- 1000 atm pressure
567   regime.  Compressibilities computed using the Hull volume, however,
568   deviate dramatically from the experimental values at low applied
569 < pressures.  The reason for this deviation is quite simple; at low
569 > pressures.  The reason for this deviation is quite simple: at low
570   applied pressures, the liquid is in equilibrium with a vapor phase,
571   and it is entirely possible for one (or a few) molecules to drift away
572   from the liquid cluster (see Fig. \ref{fig:coneOfShame}).  At low
# Line 471 | Line 575 | geometries which include large volumes of empty space.
575   geometries which include large volumes of empty space.
576  
577   \begin{figure}
578 < \includegraphics[width=\linewidth]{flytest2}
578 > \includegraphics[width=\linewidth]{coneOfShame}
579   \caption{At low pressures, the liquid is in equilibrium with the vapor
580    phase, and isolated molecules can detach from the liquid droplet.
581 <  This is expected behavior, but the reported volume of the convex
582 <  hull includes large regions of empty space.  For this reason,
581 >  This is expected behavior, but the volume of the convex hull
582 >  includes large regions of empty space. For this reason,
583    compressibilities are computed using local number densities rather
584    than hull volumes.}
585   \label{fig:coneOfShame}
586   \end{figure}
587  
588 < At higher pressures, the equilibrium favors the liquid phase, and the
589 < hull geometries are much more compact.  Because of the liquid-vapor
590 < effect on the convex hull, the regional number density approach
591 < (Eq. \ref{eq:BMN}) provides more reliable estimates of the bulk
592 < modulus.
588 > At higher pressures, the equilibrium strongly favors the liquid phase,
589 > and the hull geometries are much more compact.  Because of the
590 > liquid-vapor effect on the convex hull, the regional number density
591 > approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the
592 > compressibility.
593  
594 < We initially used the classic compressibility formula to calculate the the isothermal compressibility at each target pressure. These calculations yielded compressibility values that were dramatically higher than both previous simulations and experiment. The particular compressibility expression used requires the calculation of both a volume and pressure differential, thereby stipulating that the data from at least two simulations at different pressures must be used to calculate the isothermal compressibility at one pressure.
595 <
596 < Regardless of the difficulty in obtaining accurate hull
597 < volumes at low temperature and pressures, the Langevin Hull NPT method
598 < provides reasonable isothermal compressibility values for water
599 < through a large range of pressures.
496 <
497 < Per the fluctuation dissipation theorem \cite{Debenedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure
498 <
594 > In both the traditional compressibility formula (Eq. \ref{eq:BM}) and
595 > the number density version (Eq. \ref{eq:BMN}), multiple simulations at
596 > different pressures must be done to compute the first derivatives.  It
597 > is also possible to compute the compressibility using the fluctuation
598 > dissipation theorem using either fluctuations in the
599 > volume,\cite{Debenedetti1986}
600   \begin{equation}
601 < \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T}
601 > \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle
602 >    V \right \rangle ^{2}}{V \, k_{B} \, T},
603 > \label{eq:BMVfluct}
604   \end{equation}
605 + or, equivalently, fluctuations in the number of molecules within the
606 + fixed region,
607 + \begin{equation}
608 + \kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle
609 +    N \right \rangle ^{2}}{N \, k_{B} \, T}.
610 + \label{eq:BMNfluct}
611 + \end{equation}
612 + Thus, the compressibility of each simulation can be calculated
613 + entirely independently from other trajectories.  Compressibility
614 + calculations that rely on the hull volume will still suffer the
615 + effects of the empty space due to the vapor phase; for this reason, we
616 + recommend using the number density (Eq. \ref{eq:BMN}) or number
617 + density fluctuations (Eq. \ref{eq:BMNfluct}) for computing
618 + compressibilities.  We obtained the results in
619 + Fig. \ref{fig:compWater} using a sampling radius that was
620 + approximately 80\% of the mean distance between the center of mass of
621 + the cluster and the hull atoms.  This ratio of sampling radius to
622 + average hull radius excludes the problematic vapor phase on the
623 + outside of the cluster while including enough of the liquid phase to
624 + avoid poor statistics due to fluctuating local densities.
625  
626 < Thus, the compressibility of each simulation run can be calculated entirely independently from all other trajectories. However, the resulting compressibilities were still as much as an order of magnitude larger than the reference values. The effect was particularly pronounced at the low end of the pressure range. At ambient temperature and low pressures, there exists an equilibrium between vapor and liquid phases. Vapor molecules are naturally more diffuse around the exterior of the cluster, causing artificially large cluster volumes. Any compressibility calculation that relies on the hull volume will suffer these effects.
626 > A comparison of the oxygen-oxygen radial distribution functions for
627 > SPC/E water simulated using both the Langevin Hull and more
628 > traditional periodic boundary methods -- both at 1 atm and 300K --
629 > reveals an understructuring of water in the Langevin Hull that
630 > manifests as a slight broadening of the solvation shells.  This effect
631 > may be due to the introduction of a liquid-vapor interface in the
632 > Langevin Hull simulations (an interface which is missing in most
633 > periodic simulations of bulk water).  Vapor-phase molecules contribute
634 > a small but nearly flat portion of the radial distribution function.
635  
505 –
636   \subsection{Molecular orientation distribution at cluster boundary}
637  
638 < In order for non-periodic boundary conditions to be widely applicable, they must be constructed in such a way that they allow a finite, usually small, simulated system to replicate the properties of an infinite bulk system. Naturally, this requirement has spawned many methods for inserting boundaries into simulated systems [REF... ?]. Of particular interest to our characterization of the Langevin Hull is the orientation of water molecules included in the geometric hull. Ideally, all molecules in the cluster will have the same orientational distribution as bulk water.
638 > In order for a non-periodic boundary method to be widely applicable,
639 > it must be constructed in such a way that they allow a finite system
640 > to replicate the properties of the bulk. Early non-periodic simulation
641 > methods (e.g. hydrophobic boundary potentials) induced spurious
642 > orientational correlations deep within the simulated
643 > system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
644 > fixing and characterizing the effects of artificial boundaries
645 > including methods which fix the orientations of a set of edge
646 > molecules.\cite{Warshel1978,King1989}
647  
648 < The orientation of molecules at the edges of a simulated cluster has long been a 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, 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 70\% of the 100 molecules in the simulation. This result echoes an earlier study  which showed that an extended planar hydrophobic surface caused orientational preference at the surface which extended 7 \r{A} 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.
649 <
650 < In contrast, the Langevin Hull does not require that the orientation of molecules be fixed, nor does it utilize an explicitly hydrophobic boundary, orientational constraint or radial constraint. The number and identity of the molecules included on the convex hull are dynamic properties, thus avoiding the formation of an artificial solvent boundary layer. The hope is that the water molecules on the surface of the cluster, if left to their own devices in the absence of orientational and radial constraints, will maintain a bulk-like orientational distribution.
651 <
652 < 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.
653 <
654 < The orientation of a water molecule is described by
648 > As described above, the Langevin Hull does not require that the
649 > orientation of molecules be fixed, nor does it utilize an explicitly
650 > hydrophobic boundary, or orientational or radial constraints.
651 > Therefore, the orientational correlations of the molecules in water
652 > clusters are of particular interest in testing this method.  Ideally,
653 > the water molecules on the surfaces of the clusters will have enough
654 > mobility into and out of the center of the cluster to maintain
655 > bulk-like orientational distribution in the absence of orientational
656 > and radial constraints.  However, since the number of hydrogen bonding
657 > partners available to molecules on the exterior are limited, it is
658 > likely that there will be an effective hydrophobicity of the hull.
659  
660 + To determine the extent of these effects, we examined the
661 + orientations exhibited by SPC/E water in a cluster of 1372
662 + molecules at 300 K and at pressures ranging from 1 -- 1000 atm.  The
663 + orientational angle of a water molecule is described by
664   \begin{equation}
665   \cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
666   \end{equation}
667 + where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of
668 + mass and the cluster center of mass, and $\vec{\mu}_{i}$ is the vector
669 + bisecting the H-O-H angle of molecule {\it i}.  Bulk-like
670 + distributions will result in $\langle \cos \theta \rangle$ values
671 + close to zero.  If the hull exhibits an overabundance of
672 + externally-oriented oxygen sites, the average orientation will be
673 + negative, while dangling hydrogen sites will result in positive
674 + average orientations.
675  
676 < 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}.
677 <
676 > Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
677 > for molecules in the interior of the cluster (squares) and for
678 > molecules included in the convex hull (circles).
679   \begin{figure}
680 < \includegraphics[width=\linewidth]{g_r_theta}
681 < \caption{Definition of coordinates}
682 < \label{coords}
680 > \includegraphics[width=\linewidth]{pAngle}
681 > \caption{Distribution of $\cos{\theta}$ values for molecules on the
682 >  interior of the cluster (squares) and for those participating in the
683 >  convex hull (circles) at a variety of pressures.  The Langevin Hull
684 >  exhibits minor dewetting behavior with exposed oxygen sites on the
685 >  hull water molecules.  The orientational preference for exposed
686 >  oxygen appears to be independent of applied pressure. }
687 > \label{fig:pAngle}
688   \end{figure}
689  
690 < 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).
690 > As expected, interior molecules (those not included in the convex
691 > hull) maintain a bulk-like structure with a uniform distribution of
692 > orientations. Molecules included in the convex hull show a slight
693 > preference for values of $\cos{\theta} < 0.$ These values correspond
694 > to molecules with oxygen directed toward the exterior of the cluster,
695 > forming dangling hydrogen bond acceptor sites.
696  
697 + The orientational preference exhibited by water molecules on the hull
698 + is significantly weaker than the preference caused by an explicit
699 + hydrophobic bounding potential.  Additionally, the Langevin Hull does
700 + not require that the orientation of any molecules be fixed in order to
701 + maintain bulk-like structure, even near the cluster surface.
702 +
703 + Previous molecular dynamics simulations of SPC/E liquid / vapor
704 + interfaces using periodic boundary conditions have shown that
705 + molecules on the liquid side of interface favor a similar orientation
706 + where oxygen is directed away from the bulk.\cite{Taylor1996} These
707 + simulations had well-defined liquid and vapor phase regions
708 + equilibrium and it was observed that {\it vapor} molecules generally
709 + had one hydrogen protruding from the surface, forming a dangling
710 + hydrogen bond donor. Our water clusters do not have a true vapor
711 + region, but rather a few transient molecules that leave the liquid
712 + droplet (and which return to the droplet relatively quickly).
713 + Although we cannot obtain an orientational preference of vapor phase
714 + molecules in a Langevin Hull simulation, but we do agree with previous
715 + estimates of the orientation of {\it liquid phase} molecules at the
716 + interface.
717 +
718 + \subsection{Heterogeneous nanoparticle / water mixtures}
719 +
720 + To further test the method, we simulated gold nanoparticles ($r = 18$
721 + \AA~, 1433 atoms) solvated by explicit SPC/E water clusters (5000
722 + molecules) using a model for the gold / water interactions that has
723 + been used by Dou {\it et. al.} for investigating the separation of
724 + water films near hot metal surfaces.\cite{ISI:000167766600035} The
725 + Langevin Hull was used to sample pressures of 1, 2, 5, 10, 20, 50, 100
726 + and 200 atm, while all simulations were done at a temperature of 300
727 + K.  At these temperatures and pressures, there is no observed
728 + separation of the water film from the surface.
729 +
730 + In Fig. \ref{fig:RhoR} we show the density of water and gold as a
731 + function of the distance from the center of the nanoparticle.  Higher
732 + applied pressures appear to destroy structural correlations in the
733 + outermost monolayer of the gold nanoparticle as well as in the water
734 + at the near the metal / water interface.  Simulations at increased
735 + pressures exhibit significant overlap of the gold and water densities,
736 + indicating a less well-defined interfacial surface.
737 +
738   \begin{figure}
739 < \includegraphics[width=\linewidth]{pAngle}
740 < \caption{SPC/E water clusters: only minor dewetting at the boundary}
741 < \label{pAngle}
739 > \includegraphics[width=\linewidth]{RhoR}
740 > \caption{Density profiles of gold and water at the nanoparticle
741 >  surface. Each curve has been normalized by the average density in
742 >  the bulk-like region available to the corresponding material.
743 >  Higher applied pressures de-structure both the gold nanoparticle
744 >  surface and water at the metal/water interface.}
745 > \label{fig:RhoR}
746   \end{figure}
747  
748 < 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.
748 > At even higher pressures (500 atm and above), problems with the metal
749 > - water interaction potential became quite clear.  The model we are
750 > using appears to have been parameterized for relatively low pressures;
751 > it utilizes both shifted Morse and repulsive Morse potentials to model
752 > the Au/O and Au/H interactions, respectively.  The repulsive wall of
753 > the Morse potential does not diverge quickly enough at short distances
754 > to prevent water from diffusing into the center of the gold
755 > nanoparticles.  This behavior is likely not a realistic description of
756 > the real physics of the situation.  A better model of the gold-water
757 > adsorption behavior would require harder repulsive walls to prevent
758 > this behavior.
759  
760 < 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.
760 > \section{Discussion}
761 > \label{sec:discussion}
762  
763 < 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.
763 > The Langevin Hull samples the isobaric-isothermal ensemble for
764 > non-periodic systems by coupling the system to a bath characterized by
765 > pressure, temperature, and solvent viscosity.  This enables the
766 > simulation of heterogeneous systems composed of materials with
767 > significantly different compressibilities.  Because the boundary is
768 > dynamically determined during the simulation and the molecules
769 > interacting with the boundary can change, the method inflicts minimal
770 > perturbations on the behavior of molecules at the edges of the
771 > simulation.  Further work on this method will involve implicit
772 > electrostatics at the boundary (which is missing in the current
773 > implementation) as well as more sophisticated treatments of the
774 > surface geometry (alpha
775 > shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight
776 > Cocone\cite{Dey:2003ts}). The non-convex hull geometries are
777 > significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull
778 > ($\mathcal{O}(N \log N)$), but would enable the use of hull volumes
779 > directly in computing the compressibility of the sample.
780  
781 < \subsection{Heterogeneous nanoparticle / water mixtures}
781 > \section*{Appendix A: Computing Convex Hulls on Parallel Computers}
782  
783 + In order to use the Langevin Hull for simulations on parallel
784 + computers, one of the more difficult tasks is to compute the bounding
785 + surface, facets, and resistance tensors when the individual processors
786 + have incomplete information about the entire system's topology.  Most
787 + parallel decomposition methods assign primary responsibility for the
788 + motion of an atomic site to a single processor, and we can exploit
789 + this to efficiently compute the convex hull for the entire system.
790  
791 < \section*{Appendix A: Computing Convex Hulls on Parallel Computers}
791 > The basic idea involves splitting the point cloud into
792 > spatially-overlapping subsets and computing the convex hulls for each
793 > of the subsets.  The points on the convex hull of the entire system
794 > are all present on at least one of the subset hulls. The algorithm
795 > works as follows:
796 > \begin{enumerate}
797 > \item Each processor computes the convex hull for its own atomic sites
798 >  (left panel in Fig. \ref{fig:parallel}).
799 > \item The Hull vertices from each processor are communicated to all of
800 >  the processors, and each processor assembles a complete list of hull
801 >  sites (this is much smaller than the original number of points in
802 >  the point cloud).
803 > \item Each processor computes the global convex hull (right panel in
804 >  Fig. \ref{fig:parallel}) using only those points that are the union
805 >  of sites gathered from all of the subset hulls.  Delaunay
806 >  triangulation is then done to obtain the facets of the global hull.
807 > \end{enumerate}
808  
809 + \begin{figure}
810 + \includegraphics[width=\linewidth]{parallel}
811 + \caption{When the sites are distributed among many nodes for parallel
812 +  computation, the processors first compute the convex hulls for their
813 +  own sites (dashed lines in left panel). The positions of the sites
814 +  that make up the subset hulls are then communicated to all
815 +  processors (middle panel).  The convex hull of the system (solid line in
816 +  right panel) is the convex hull of the points on the union of the subset
817 +  hulls.}
818 + \label{fig:parallel}
819 + \end{figure}
820 +
821 + The individual hull operations scale with
822 + $\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total
823 + number of sites, and $p$ is the number of processors.  These local
824 + hull operations create a set of $p$ hulls, each with approximately
825 + $\frac{n}{3pr}$ sites for a cluster of radius $r$. The worst-case
826 + communication cost for using a ``gather'' operation to distribute this
827 + information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n
828 +  \beta (p-1)}{3 r p^2})$, while the final computation of the system
829 + hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$.
830 +
831 + For a large number of atoms on a moderately parallel machine, the
832 + total costs are dominated by the computations of the individual hulls,
833 + and communication of these hulls to create the Langevin Hull sees roughly
834 + linear speed-up with increasing processor counts.
835 +
836   \section*{Acknowledgments}
837   Support for this project was provided by the
838   National Science Foundation under grant CHE-0848243. Computational
839   time was provided by the Center for Research Computing (CRC) at the
840   University of Notre Dame.  
841  
842 + Molecular graphics images were produced using the UCSF Chimera package from
843 + the Resource for Biocomputing, Visualization, and Informatics at the
844 + University of California, San Francisco (supported by NIH P41 RR001081).
845   \newpage
846  
847   \bibliography{langevinHull}

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines