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 3652 by gezelter, Mon Oct 18 18:27:24 2010 UTC vs.
Revision 3698 by kstocke1, Wed Nov 17 19:22:20 2010 UTC

# Line 18 | Line 18
18   \setlength{\belowcaptionskip}{30 pt}
19  
20   \bibpunct{[}{]}{,}{s}{}{;}
21 < \bibliographystyle{aip}
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.
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.
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
# 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.
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 < As long as the material in the simulation box is essentially a bulk
91 < liquid which has a relatively uniform compressibility, the standard
92 < approach provides an excellent way of adjusting the volume of the
93 < system and applying pressure directly via the interactions between
94 < atomic sites.  
90 <
91 < The problem with these approaches becomes apparent when the material
92 < being simulated is an inhomogeneous mixture in which portions of the
93 < simulation box are incompressible relative to other portions.
94 < Examples include simulations of metallic nanoparticles in liquid
95 < 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 collisions 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 < Additionally, one may often wish to simulate explicitly non-periodic
114 < systems, and the constraint that a periodic box must be used to
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 < Explicitly non-periodic systems
124 > {\it Yotal} 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 < Elastic Bag
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 < Spherical Boundary approaches
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 < \section{Methodology}
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 < We have developed a new method which uses a constant pressure and
172 < temperature bath.  This bath interacts only with the objects that are
173 < currently at the edge of the system.  Since the edge is determined
174 < dynamically as the simulation progresses, no {\it a priori} geometry
175 < is defined.  The pressure and temperature bath interacts {\it
176 <  directly} with the atoms on the edge and not with atoms interior to
177 < the simulation.  This means that there are no affine transforms
178 < required.  There are also no fictitious particles or bounding
179 < potentials used in this approach.
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 < The basics of the method are as follows. The simulation starts as a
182 < collection of atomic locations in three dimensions (a point cloud).
183 < Delaunay triangulation is used to find all facets between coplanar
184 < neighbors.  In highly symmetric point clouds, facets can contain many
185 < atoms, but in all but the most symmetric of cases one might experience
186 < in a molecular dynamics simulation, the facets are simple triangles in
187 < 3-space that contain exactly three atoms.  
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 < The convex hull is the set of facets that have {\it no concave
195 <  corners} at an atomic site.  This eliminates all facets on the
196 < interior of the point cloud, leaving only those exposed to the
197 < bath. Sites on the convex hull are dynamic. As molecules re-enter the
198 < cluster, all interactions between atoms on that molecule and the
199 < external bath are removed.
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 < For atomic sites in the interior of the point cloud, the equations of
211 < motion are simple Newtonian dynamics,
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}
# Line 158 | 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 176 | Line 285 | coupling depends on the solvent temperature, friction
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, friction and the size and
289 < shape of each facet. The thermal interactions are expressed as a
290 < typical Langevin description of the forces,
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, $P$ is the external pressure, $A_f$ and $\hat{n}_f$ are the area
298 < and normal vectors for facet $f$, respectively.  ${\mathbf v}_f(t)$ is
299 < the velocity of the facet,
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 a ($3 \times 3$) hydrodynamic tensor that depends on
304 < the geometry and surface area of facet $f$ and the viscosity of the
305 < fluid (See Appendix A).  The hydrodynamic tensor is related to the
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}
# Line 203 | Line 312 | Once the hydrodynamic tensor is known for a given face
312   \label{eq:randomForce}
313   \end{eqnarray}
314  
315 < Once the hydrodynamic tensor is known for a given facet (see Appendix
316 < A) obtaining a stochastic vector that has the properties in
317 < Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a
318 < one-time Cholesky decomposition to obtain the square root matrix of
210 < the resistance tensor,
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}
# Line 224 | Line 332 | We have implemented this method by extending the Lange
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 < We have implemented this method by extending the Langevin dynamics
336 < integrator in our group code, OpenMD.\cite{Meineke2005,openmd}  
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 < \section{Tests \& Applications}
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  
232 \subsection{Bulk modulus of gold nanoparticles}
233
358   \begin{figure}
359 < \includegraphics[width=\linewidth]{pressure_tb}
360 < \caption{Pressure response is rapid (18 \AA gold nanoparticle), target
361 < pressure = 4 GPa}
362 < \label{pressureResponse}
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 < \begin{figure}
370 < \includegraphics[width=\linewidth]{temperature_tb}
371 < \caption{Temperature equilibration depends on surface area and bath
244 <  viscosity.  Target Temperature = 300K}
245 < \label{temperatureResponse}
246 < \end{figure}
247 <
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 < \kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial
250 <    P}\right)
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 < \begin{figure}
381 < \includegraphics[width=\linewidth]{compress_tb}
382 < \caption{Isothermal Compressibility (18 \AA gold nanoparticle)}
383 < \label{temperatureResponse}
384 < \end{figure}
385 <
386 < \subsection{Compressibility of SPC/E water clusters}
387 <
388 < Both NVT \cite{Glattli2002} and NPT \cite{Motakabbir1990, Pi2009} molecular dynamics simulations of SPC/E water have yielded values for the isothermal compressibility of water that agree well with experiment \cite{Fine1973}. The results of three different methods for computing the isothermal compressibility from Langevin Hull simulations for pressures between 1 and 6500 atm are shown in Fig. 5 along with compressibility values obtained from both other SPC/E simulations and experiment. Compressibility values from all references are for applied pressures within the range 1 - 1000 atm.
389 <
390 < \begin{figure}
391 < \includegraphics[width=\linewidth]{new_isothermal}
392 < \caption{Compressibility of SPC/E water}
393 < \label{compWater}
394 < \end{figure}
395 <
396 < We initially used the classic compressibility formula
380 > We have implemented this method by extending the Langevin dynamics
381 > integrator in our code, OpenMD.\cite{Meineke2005,openmd}  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{Qhull} 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 < \begin{equation}
410 < \kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T}
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 < 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.
423 <
424 < 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
425 <
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{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, 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
436 + found within this region throughout a given simulation. The geometry
437 + and size of the region is arbitrary, and any bulk-like portion of the
438 + cluster can be used to compute the compressibility.
439  
440 < 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.
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 < In order to calculate the isothermal compressibility without being hindered by hull volume issues, we adapted the classic compressibility formula so that the compressibility could be calculated using information about the local density instead of the volume of the convex hull. We calculated the $g_{OO}(r)$ for a 1 nanosecond simulation of a cluster of 1372 SPC/E water molecules and spherically integrated the function over the bounds 0 to $r'$. In all cases, the value of $r'$ was 17.26216 $\AA$. The value of the total integral between these bounds is essentially the number (N) of molecules within volume $\frac{4}{3}\pi r'^{3}$ at a given pressure. To yield an actual molecule count, N must be scaled by an ideal density. However, even in the absence of an ideal density, we can use the relationship $\rho = \frac{N}{V}$ to rewrite the isothermal compressibility formula as
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 < \kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T}
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,QSC}
493  
494 < Isothermal compressibility values calculated using this modified expression are in good agreement with the reference values throughout the 1 - 1000 atm pressure regime. Regardless of the difficulty in obtaining accurate hull volumes at low temperature and pressures, the Langevin Hull NPT method provides reasonable isothermal compressibility values for water through a large range of pressures.
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{QSC} Using the same force field, we have performed
498 > a series of 1 ns simulations on 40 \AA~ radius
499 > nanoparticles under the Langevin Hull at a variety of applied
500 > pressures ranging from 0 -- 10 GPa.  We obtain a value of 177.55 GPa
501 > for the bulk modulus of gold using this technique, in close agreement
502 > with both previous simulations and the experimental bulk modulus of
503 > gold.
504  
505 < \subsection{Molecular orientation distribution at cluster boundary}
505 > \begin{figure}
506 > \includegraphics[width=\linewidth]{stacked}
507 > \caption{The response of the internal pressure and temperature of gold
508 >  nanoparticles when first placed in the Langevin Hull
509 >  ($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting
510 >  from initial conditions that were far from the bath pressure and
511 >  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).}
512 > \label{fig:pressureResponse}
513 > \end{figure}
514  
515 < 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.
515 > We note that the Langevin Hull produces rapidly-converging behavior
516 > for structures that are started far from equilibrium.  In
517 > Fig. \ref{fig:pressureResponse} we show how the pressure and
518 > temperature respond to the Langevin Hull for nanoparticles that were
519 > initialized far from the target pressure and temperature.  As
520 > expected, the rate at which thermal equilibrium is achieved depends on
521 > the total surface area of the cluter exposed to the bath as well as
522 > the bath viscosity.  Pressure that is applied suddenly to a cluster
523 > can excite breathing vibrations, but these rapidly damp out (on time
524 > scales of 30 -- 50 ps).
525  
526 < 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.
526 > \subsection{Compressibility of SPC/E water clusters}
527  
528 < 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.
528 > Prior molecular dynamics simulations on SPC/E water (both in
529 > NVT~\cite{Glattli2002} and NPT~\cite{Motakabbir1990, Pi2009}
530 > ensembles) have yielded values for the isothermal compressibility that
531 > agree well with experiment.\cite{Fine1973} The results of two
532 > different approaches for computing the isothermal compressibility from
533 > Langevin Hull simulations for pressures between 1 and 6500 atm are
534 > shown in Fig. \ref{fig:compWater} along with compressibility values
535 > obtained from both other SPC/E simulations and experiment.
536  
537 < 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.
537 > \begin{figure}
538 > \includegraphics[width=\linewidth]{new_isothermalN}
539 > \caption{Compressibility of SPC/E water}
540 > \label{fig:compWater}
541 > \end{figure}
542  
543 < The orientation of a water molecule is described by
543 > Isothermal compressibility values calculated using the number density
544 > (Eq. \ref{eq:BMN}) expression are in good agreement with experimental
545 > and previous simulation work throughout the 1 -- 1000 atm pressure
546 > regime.  Compressibilities computed using the Hull volume, however,
547 > deviate dramatically from the experimental values at low applied
548 > pressures.  The reason for this deviation is quite simple; at low
549 > applied pressures, the liquid is in equilibrium with a vapor phase,
550 > and it is entirely possible for one (or a few) molecules to drift away
551 > from the liquid cluster (see Fig. \ref{fig:coneOfShame}).  At low
552 > pressures, the restoring forces on the facets are very gentle, and
553 > this means that the hulls often take on relatively distorted
554 > geometries which include large volumes of empty space.
555  
556 + \begin{figure}
557 + \includegraphics[width=\linewidth]{coneOfShame}
558 + \caption{At low pressures, the liquid is in equilibrium with the vapor
559 +  phase, and isolated molecules can detach from the liquid droplet.
560 +  This is expected behavior, but the volume of the convex hull
561 +  includes large regions of empty space. For this reason,
562 +  compressibilities are computed using local number densities rather
563 +  than hull volumes.}
564 + \label{fig:coneOfShame}
565 + \end{figure}
566 +
567 + At higher pressures, the equilibrium strongly favors the liquid phase,
568 + and the hull geometries are much more compact.  Because of the
569 + liquid-vapor effect on the convex hull, the regional number density
570 + approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the
571 + compressibility.
572 +
573 + In both the traditional compressibility formula (Eq. \ref{eq:BM}) and
574 + the number density version (Eq. \ref{eq:BMN}), multiple simulations at
575 + different pressures must be done to compute the first derivatives.  It
576 + is also possible to compute the compressibility using the fluctuation
577 + dissipation theorem using either fluctuations in the
578 + volume,\cite{Debenedetti1986},
579 + \begin{equation}
580 + \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle
581 +    V \right \rangle ^{2}}{V \, k_{B} \, T},
582 + \label{eq:BMVfluct}
583 + \end{equation}
584 + or, equivalently, fluctuations in the number of molecules within the
585 + fixed region,
586 + \begin{equation}
587 + \kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle
588 +    N \right \rangle ^{2}}{N \, k_{B} \, T}.
589 + \label{eq:BMNfluct}
590 + \end{equation}
591 + Thus, the compressibility of each simulation can be calculated
592 + entirely independently from other trajectories.  Compressibility
593 + calculations that rely on the hull volume will still suffer the
594 + effects of the empty space due to the vapor phase; for this reason, we
595 + recommend using the number density (Eq. \ref{eq:BMN}) or number
596 + density fluctuations (Eq. \ref{eq:BMNfluct}) for computing
597 + compressibilities.
598 +
599 + \subsection{Molecular orientation distribution at cluster boundary}
600 +
601 + In order for a non-periodic boundary method to be widely applicable,
602 + it must be constructed in such a way that they allow a finite system
603 + to replicate the properties of the bulk. Early non-periodic simulation
604 + methods (e.g. hydrophobic boundary potentials) induced spurious
605 + orientational correlations deep within the simulated
606 + system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
607 + fixing and characterizing the effects of artifical boundaries
608 + including methods which fix the orientations of a set of edge
609 + molecules.\cite{Warshel1978,King1989}
610 +
611 + As described above, the Langevin Hull does not require that the
612 + orientation of molecules be fixed, nor does it utilize an explicitly
613 + hydrophobic boundary, or orientational or radial constraints.
614 + Therefore, the orientational correlations of the molecules in water
615 + clusters are of particular interest in testing this method.  Ideally,
616 + the water molecules on the surfaces of the clusterss will have enough
617 + mobility into and out of the center of the cluster to maintain
618 + bulk-like orientational distribution in the absence of orientational
619 + and radial constraints.  However, since the number of hydrogen bonding
620 + partners available to molecules on the exterior are limited, it is
621 + likely that there will be an effective hydrophobicity of the hull.
622 +
623 + To determine the extent of these effects, we examined the
624 + orientations exhibited by SPC/E water in a cluster of 1372
625 + molecules at 300 K and at pressures ranging from 1 -- 1000 atm.  The
626 + orientational angle of a water molecule is described by
627   \begin{equation}
628   \cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
629   \end{equation}
630 + where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of
631 + mass and the cluster center of mass, and $\vec{\mu}_{i}$ is the vector
632 + bisecting the H-O-H angle of molecule {\it i}.  Bulk-like
633 + distributions will result in $\langle \cos \theta \rangle$ values
634 + close to zero.  If the hull exhibits an overabundance of
635 + externally-oriented oxygen sites, the average orientation will be
636 + negative, while dangling hydrogen sites will result in positive
637 + average orientations.
638  
639 < 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}.
640 <
639 > Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
640 > for molecules in the interior of the cluster (squares) and for
641 > molecules included in the convex hull (circles).
642   \begin{figure}
643 < \includegraphics[width=\linewidth]{g_r_theta}
644 < \caption{Definition of coordinates}
645 < \label{coords}
643 > \includegraphics[width=\linewidth]{pAngle}
644 > \caption{Distribution of $\cos{\theta}$ values for molecules on the
645 >  interior of the cluster (squares) and for those participating in the
646 >  convex hull (circles) at a variety of pressures.  The Langevin Hull
647 >  exhibits minor dewetting behavior with exposed oxygen sites on the
648 >  hull water molecules.  The orientational preference for exposed
649 >  oxygen appears to be independent of applied pressure. }
650 > \label{fig:pAngle}
651   \end{figure}
652  
653 < 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).
653 > As expected, interior molecules (those not included in the convex
654 > hull) maintain a bulk-like structure with a uniform distribution of
655 > orientations. Molecules included in the convex hull show a slight
656 > preference for values of $\cos{\theta} < 0.$ These values correspond
657 > to molecules with oxygen directed toward the exterior of the cluster,
658 > forming a dangling hydrogen bond acceptor site.
659  
660 + The orientational preference exhibited by liquid phase hull molecules in the Langevin Hull is significantly weaker than the preference caused by an explicit hydrophobic bounding potential.  Additionally, the Langevin Hull does not require that the orientation of any molecules be fixed in order to maintain bulk-like structure, even at the cluster surface.
661 +
662 + Previous molecular dynamics simulations
663 + of SPC/E water using periodic boundary conditions have shown that molecules on the liquid side of the liquid/vapor interface favor a similar orientation where oxygen is directed away from the bulk.\cite{Taylor1996} These simulations had both a liquid phase and a well-defined vapor phase in equilibrium and showed that vapor molecules generally had one hydrogen protruding from the surface, forming a dangling hydrogen bond donor. Our water cluster simulations do not have a true lasting vapor phase, but rather a few transient molecules that leave the liquid droplet. Thus while we are unable to comment on the orientational preference of vapor phase molecules in a Langevin Hull simulation, we achieve good agreement for the orientation of liquid phase molecules at the interface.
664 +
665 + \subsection{Heterogeneous nanoparticle / water mixtures}
666 +
667 + To further test the method, we simulated gold nanopartices ($r = 18$
668 + \AA) solvated by explicit SPC/E water clusters using the Langevin
669 + Hull.  This was done at pressures of 1, 2, 5, 10, 20, 50, 100 and 200 atm
670 + in order to observe the effects of pressure on the ordering of water
671 + ordering at the surface.  In Fig. \ref{fig:RhoR} we show the density
672 + of water adjacent to the surface as a function of pressure, as well as
673 + the orientational ordering of water at the surface of the
674 + nanoparticle.
675 +
676   \begin{figure}
677 < \includegraphics[width=\linewidth]{pAngle}
678 < \caption{SPC/E water clusters: only minor dewetting at the boundary}
679 < \label{pAngle}
677 >
678 > \caption{Higher applied pressures de-structure both the gold nanoparticle and water at the metal/water interface.}
679 > \label{fig:RhoR}
680   \end{figure}
681  
682 < 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.
682 > At higher pressures, problems with the gold - water interaction
683 > potential became apparent.  The model we are using (due to Spohr) was
684 > intended for relatively low pressures; it utilizes both shifted Morse
685 > and repulsive Morse potentials to model the Au/O and Au/H
686 > interactions, respectively.  The repulsive wall of the Morse potential
687 > does not diverge quickly enough at short distances to prevent water
688 > from diffusing into the center of the gold nanoparticles.  This
689 > behavior is likely not a realistic description of the real physics of
690 > the situation.  A better model of the gold-water adsorption behavior
691 > appears to require harder repulsive walls to prevent this behavior.
692  
693 < 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.
693 > \section{Discussion}
694 > \label{sec:discussion}
695  
696 < 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.
696 > The Langevin Hull samples the isobaric-isothermal ensemble for
697 > non-periodic systems by coupling the system to a bath characterized by
698 > pressure, temperature, and solvent viscosity.  This enables the
699 > simulation of heterogeneous systems composed of materials with
700 > significantly different compressibilities.  Because the boundary is
701 > dynamically determined during the simulation and the molecules
702 > interacting with the boundary can change, the method inflicts minimal
703 > perturbations on the behavior of molecules at the edges of the
704 > simulation.  Further work on this method will involve implicit
705 > electrostatics at the boundary (which is missing in the current
706 > implementation) as well as more sophisticated treatments of the
707 > surface geometry (alpha
708 > shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight
709 > Cocone\cite{Dey:2003ts}). The non-convex hull geometries are
710 > significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull
711 > ($\mathcal{O}(N \log N)$), but would enable the use of hull volumes
712 > directly in computing the compressibility of the sample.
713  
714 + \section*{Appendix A: Computing Convex Hulls on Parallel Computers}
715  
716 < \subsection{Heterogeneous nanoparticle / water mixtures}
716 > In order to use the Langevin Hull for simulations on parallel
717 > computers, one of the more difficult tasks is to compute the bounding
718 > surface, facets, and resistance tensors when the individual processors
719 > have incomplete information about the entire system's topology.  Most
720 > parallel decomposition methods assign primary responsibility for the
721 > motion of an atomic site to a single processor, and we can exploit
722 > this to efficiently compute the convex hull for the entire system.
723  
724 + The basic idea involves splitting the point cloud into
725 + spatially-overlapping subsets and computing the convex hulls for each
726 + of the subsets.  The points on the convex hull of the entire system
727 + are all present on at least one of the subset hulls. The algorithm
728 + works as follows:
729 + \begin{enumerate}
730 + \item Each processor computes the convex hull for its own atomic sites
731 +  (left panel in Fig. \ref{fig:parallel}).
732 + \item The Hull vertices from each processor are communicated to all of
733 +  the processors, and each processor assembles a complete list of hull
734 +  sites (this is much smaller than the original number of points in
735 +  the point cloud).
736 + \item Each processor computes the global convex hull (right panel in
737 +  Fig. \ref{fig:parallel}) using only those points that are the union
738 +  of sites gathered from all of the subset hulls.  Delaunay
739 +  triangulation is then done to obtain the facets of the global hull.
740 + \end{enumerate}
741  
335 \section{Appendix A: Hydrodynamic tensor for triangular facets}
336
742   \begin{figure}
743 < \includegraphics[width=\linewidth]{hydro}
744 < \caption{Hydro}
745 < \label{hydro}
743 > \includegraphics[width=\linewidth]{parallel}
744 > \caption{When the sites are distributed among many nodes for parallel
745 >  computation, the processors first compute the convex hulls for their
746 >  own sites (dashed lines in left panel). The positions of the sites
747 >  that make up the subset hulls are then communicated to all
748 >  processors (middle panel).  The convex hull of the system (solid line in
749 >  right panel) is the convex hull of the points on the union of the subset
750 >  hulls.}
751 > \label{fig:parallel}
752   \end{figure}
753  
754 < \begin{equation}
755 < \Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}
756 < \end{equation}
754 > The individual hull operations scale with
755 > $\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total
756 > number of sites, and $p$ is the number of processors.  These local
757 > hull operations create a set of $p$ hulls, each with approximately
758 > $\frac{n}{3pr}$ sites for a cluster of radius $r$. The worst-case
759 > communication cost for using a ``gather'' operation to distribute this
760 > information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n
761 >  \beta (p-1)}{3 r p^2})$, while the final computation of the system
762 > hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$.
763  
764 < \begin{equation}
765 < T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I +
766 <  \frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right)
767 < \end{equation}
764 > For a large number of atoms on a moderately parallel machine, the
765 > total costs are dominated by the computations of the individual hulls,
766 > and communication of these hulls to create the Langevin Hull sees roughly
767 > linear speed-up with increasing processor counts.
768  
769 < \section{Appendix B: Computing Convex Hulls on Parallel Computers}
353 <
354 < \section{Acknowledgments}
769 > \section*{Acknowledgments}
770   Support for this project was provided by the
771   National Science Foundation under grant CHE-0848243. Computational
772   time was provided by the Center for Research Computing (CRC) at the
773   University of Notre Dame.  
774  
775 + Molecular graphics images were produced using the UCSF Chimera package from
776 + the Resource for Biocomputing, Visualization, and Informatics at the
777 + University of California, San Francisco (supported by NIH P41 RR001081).
778   \newpage
779  
780   \bibliography{langevinHull}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines