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 3649 by kstocke1, Thu Sep 23 21:42:37 2010 UTC vs.
Revision 3660 by gezelter, Thu Oct 21 15:27:57 2010 UTC

# Line 41 | Line 41 | Notre Dame, Indiana 46556}
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'',
45 <  performs better than traditional affine transform methods for
46 <  systems containing heterogeneous mixtures of materials with
47 <  different compressibilities. It does not suffer from the edge
48 <  effects of boundary potential methods, and allows realistic
49 <  treatment of both external pressure and thermal conductivity to an
50 <  implicit solvents.  We apply this method to several different
51 <  systems including bare nano-particles, nano-particles in explicit
52 <  solvent, as well as clusters of liquid water and ice. The predicted
53 <  mechanical and thermal properties of these systems are in good
54 <  agreement with experimental data.
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.
55   \end{abstract}
56  
57   \newpage
# Line 75 | Line 75 | Melchionna modification\cite{melchionna93} to 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
77   particles). The most common constant pressure methods, including the
78 < Melchionna modification\cite{melchionna93} to the
79 < Nos\'e-Hoover-Andersen equations of motion, the Berendsen pressure
80 < bath, and the Langevin Piston, all utilize coordinate transformation
81 < to adjust the box volume.
78 > Melchionna modification\cite{Melchionna1993} to the
79 > Nos\'e-Hoover-Andersen equations of
80 > motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen
81 > pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin
82 > Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize coordinate
83 > transformation to adjust the box volume.  As long as the material in
84 > the simulation box is essentially a bulk-like liquid which has a
85 > relatively uniform compressibility, the standard affine transform
86 > approach provides an excellent way of adjusting the volume of the
87 > system and applying pressure directly via the interactions between
88 > 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
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.
100 +
101   \begin{figure}
102   \includegraphics[width=\linewidth]{AffineScale2}
103   \caption{Affine Scaling constant pressure methods use box-length
# Line 92 | Line 110 | to adjust the box volume.
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
117 + higher than desirable, or unreasonable system sizes to avoid this
118 + effect.  For example, calculations using typical hydration shells
119 + solvating a protein under periodic boundary conditions are quite
120 + expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL
121 + SIMULATIONS]
122  
123 < Heterogeneous mixtures of materials with different compressibilities?
123 > There have been a number of other approaches to explicit
124 > non-periodicity that focus on constant or nearly-constant {\it volume}
125 > conditions while maintaining bulk-like behavior.  Berkowitz and
126 > McCammon introduced a stochastic (Langevin) boundary layer inside a
127 > region of fixed molecules which effectively enforces constant
128 > 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]
134  
135 < Explicitly non-periodic systems
135 > The electrostatic and dispersive behavior near the boundary has long
136 > been a cause for concern.  King and Warshel introduced a surface
137 > constrained all-atom solvent (SCAAS) which included polarization
138 > effects of a fixed spherical boundary to mimic bulk-like behavior
139 > without periodic boundaries.\cite{king:3647} In the SCAAS model, a
140 > layer of fixed solvent molecules surrounds the solute and any explicit
141 > solvent, and this in turn is surrounded by a continuum dielectric.
142 > MORE HERE.  WHAT DID THEY FIND?
143  
144 < Elastic Bag
144 > Beglov and Roux developed a boundary model in which the hard sphere
145 > boundary has a radius that varies with the instantaneous configuration
146 > of the solute (and solvent) molecules.\cite{beglov:9050} This model
147 > contains a clear pressure and surface tension contribution to the free
148 > energy which XXX.
149  
150 < Spherical Boundary approaches
150 > Restraining {\it potentials} introduce repulsive potentials at the
151 > surface of a sphere or other geometry.  The solute and any explicit
152 > solvent are therefore restrained inside this potential.  Often the
153 > potentials include a weak short-range attraction to maintain the
154 > correct density at the boundary.  Beglov and Roux have also introduced
155 > a restraining boundary potential which relaxes dynamically depending
156 > on the solute geometry and the force the explicit system exerts on the
157 > shell.\cite{Beglov:1995fk}
158  
159 < \section{Methodology}
159 > Recently, Krilov {\it et al.} introduced a flexible boundary model
160 > that uses a Lennard-Jones potential between the solvent molecules and
161 > a boundary which is determined dynamically from the position of the
162 > nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This approach allows
163 > the confining potential to prevent solvent molecules from migrating
164 > too far from the solute surface, while providing a weak attractive
165 > force pulling the solvent molecules towards a fictitious bulk solvent.
166 > Although this approach is appealing and has physical motivation,
167 > nanoparticles do not deform far from their original geometries even at
168 > temperatures which vaporize the nearby solvent. For the systems like
169 > the one described, the flexible boundary model will be nearly
170 > identical to a fixed-volume restraining potential.
171  
172 < A new method which uses a constant pressure and temperature bath that
173 < interacts with the objects that are currently at the edge of the
174 < system.
172 > The approach of Kohanoff, Caro, and Finnis is the most promising of
173 > the methods for introducing both constant pressure and temperature
174 > into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru}
175 > This method is based on standard Langevin dynamics, but the Brownian
176 > or random forces are allowed to act only on peripheral atoms and exert
177 > force in a direction that is inward-facing relative to the facets of a
178 > closed bounding surface.  The statistical distribution of the random
179 > forces are uniquely tied to the pressure in the external reservoir, so
180 > the method can be shown to sample the isobaric-isothermal ensemble.
181 > Kohanoff {\it et al.} used a Delaunay tessellation to generate a
182 > bounding surface surrounding the outermost atoms in the simulated
183 > system.  This is not the only possible triangulated outer surface, but
184 > guarantees that all of the random forces point inward towards the
185 > cluster.
186  
187 < Novel features: No a priori geometry is defined, No affine transforms,
188 < No fictitious particles, No bounding potentials.
187 > In the following sections, we extend and generalize the approach of
188 > Kohanoff, Caro, and Finnis. The new method, which we are calling the
189 > ``Langevin Hull'' applies the external pressure, Langevin drag, and
190 > random forces on the facets of the {\it hull itself} instead of the
191 > atomic sites comprising the vertices of the hull.  This allows us to
192 > decouple the external pressure contribution from the drag and random
193 > force.  Section \ref{sec:meth}
194  
195 < Simulation starts as a collection of atomic locations in 3D (a point
196 < cloud).
195 > \section{Methodology}
196 > \label{sec:meth}
197  
198 < Delaunay triangulation finds all facets between coplanar neighbors.
198 > We have developed a new method which uses an external bath at a fixed
199 > constant pressure ($P$) and temperature ($T$).  This bath interacts
200 > only with the objects on the exterior hull of the system.  Defining
201 > the hull of the simulation is done in a manner similar to the approach
202 > of Kohanoff, Caro and Finnis.\cite{Kohanoff:2005qm} That is, any
203 > instantaneous configuration of the atoms in the system is considered
204 > as a point cloud in three dimensional space.  Delaunay triangulation
205 > is used to find all facets between coplanar neighbors.\cite{DELAUNAY}
206 > In highly symmetric point clouds, facets can contain many atoms, but
207 > in all but the most symmetric of cases the facets are simple triangles
208 > in 3-space that contain exactly three atoms.
209  
210 < The Convex Hull is the set of facets that have no concave corners at a
211 < vertex.
210 > The convex hull is the set of facets that have {\it no concave
211 >  corners} at an atomic site.  This eliminates all facets on the
212 > interior of the point cloud, leaving only those exposed to the
213 > bath. Sites on the convex hull are dynamic. As molecules re-enter the
214 > cluster, all interactions between atoms on that molecule and the
215 > external bath are removed.  Since the edge is determined dynamically
216 > as the simulation progresses, no {\it a priori} geometry is
217 > defined. The pressure and temperature bath interacts {\it directly}
218 > with the atoms on the edge and not with atoms interior to the
219 > simulation.
220  
221 < Molecules on the convex hull are dynamic. As they re-enter the
222 < cluster, all interactions with the external bath are removed.The
123 < external bath applies pressure to the facets of the convex hull in
124 < direct proportion to the area of the facet. Thermal coupling depends on
125 < the solvent temperature, friction and the size and shape of each
126 < facet.
127 <
221 > Atomic sites in the interior of the point cloud move under standard
222 > Newtonian dynamics,
223   \begin{equation}
224 < m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U
224 > m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U,
225 > \label{eq:Newton}
226   \end{equation}
227 <
227 > where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the
228 > instantaneous velocity of site $i$ at time $t$, and $U$ is the total
229 > potential energy.  For atoms on the exterior of the cluster
230 > (i.e. those that occupy one of the vertices of the convex hull), the
231 > equation of motion is modified with an external force, ${\mathbf
232 >  F}_i^{\mathrm ext}$,
233   \begin{equation}
234 < m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}
234 > m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}.
235   \end{equation}
236  
237 + The external bath interacts directly with the facets of the convex
238 + hull. Since each vertex (or atom) provides one corner of a triangular
239 + facet, the force on the facets are divided equally to each vertex.
240 + However, each vertex can participate in multiple facets, so the resultant
241 + force is a sum over all facets $f$ containing vertex $i$:
242   \begin{equation}
243   {\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\
244      } f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\  {\mathbf
245    F}_f^{\mathrm ext}
246   \end{equation}
247  
248 + The external pressure bath applies a force to the facets of the convex
249 + hull in direct proportion to the area of the facet, while the thermal
250 + coupling depends on the solvent temperature, viscosity and the size
251 + and shape of each facet. The thermal interactions are expressed as a
252 + standard Langevin description of the forces,
253   \begin{equation}
254   \begin{array}{rclclcl}
255   {\mathbf F}_f^{\text{ext}} & = &  \text{external pressure} & + & \text{drag force} & + & \text{random force} \\
256   & = &  -\hat{n}_f P A_f  & - & \Xi_f(t) {\mathbf v}_f(t)  & + & {\mathbf R}_f(t)
257   \end{array}
258   \end{equation}
259 <
260 < \begin{eqnarray}
261 < A_f & = & \text{area of facet}\ f \\
262 < \hat{n}_f & = & \text{facet normal} \\
263 < P & = & \text{external pressure}
264 < \end{eqnarray}
265 <
266 < \begin{eqnarray}
267 < {\mathbf v}_f(t) & = & \text{velocity of facet} \\
268 < & = & \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i \\
269 < \Xi_f(t) & = & \text{is a hydrodynamic tensor that depends} \\
159 < & & \text{on the geometry and surface area of} \\
160 < & & \text{facet} \ f\ \text{and the viscosity of the fluid.}
161 < \end{eqnarray}
162 <
259 > Here, $A_f$ and $\hat{n}_f$ are the area and normal vectors for facet
260 > $f$, respectively.  ${\mathbf v}_f(t)$ is the velocity of the facet
261 > centroid,
262 > \begin{equation}
263 > {\mathbf v}_f(t) =  \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i,
264 > \end{equation}
265 > and $\Xi_f(t)$ is an approximate ($3 \times 3$) resistance tensor that
266 > depends on the geometry and surface area of facet $f$ and the
267 > viscosity of the fluid.  The resistance tensor is related to the
268 > fluctuations of the random force, $\mathbf{R}(t)$, by the
269 > fluctuation-dissipation theorem,
270   \begin{eqnarray}
271   \left< {\mathbf R}_f(t) \right> & = & 0 \\
272   \left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\
273 < \Xi_f(t)\delta(t-t^\prime)
273 > \Xi_f(t)\delta(t-t^\prime).
274 > \label{eq:randomForce}
275   \end{eqnarray}
276  
277 < Implemented in OpenMD.\cite{Meineke:2005gd,openmd}
277 > Once the resistance tensor is known for a given facet a stochastic
278 > vector that has the properties in Eq. (\ref{eq:randomForce}) can be
279 > done efficiently by carrying out a Cholesky decomposition to obtain
280 > the square root matrix of the resistance tensor,
281 > \begin{equation}
282 > \Xi_f = {\bf S} {\bf S}^{T},
283 > \label{eq:Cholesky}
284 > \end{equation}
285 > where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
286 > vector with the statistics required for the random force can then be
287 > obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which
288 > has elements chosen from a Gaussian distribution, such that:
289 > \begin{equation}
290 > \langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
291 > {\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
292 > \end{equation}
293 > where $\delta t$ is the timestep in use during the simulation. The
294 > random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to
295 > have the correct properties required by Eq. (\ref{eq:randomForce}).
296  
297 + Our treatment of the resistance tensor is approximate.  $\Xi$ for a
298 + rigid triangular plate would normally be treated as a $6 \times 6$
299 + tensor that includes translational and rotational drag as well as
300 + translational-rotational coupling. The computation of resistance
301 + tensors for rigid bodies has been detailed
302 + elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun2008}
303 + but the standard approach involving bead approximations would be
304 + prohibitively expensive if it were recomputed at each step in a
305 + molecular dynamics simulation.
306 +
307 + We are utilizing an approximate resistance tensor obtained by first
308 + constructing the Oseen tensor for the interaction of the centroid of
309 + the facet ($f$) with each of the subfacets $j$,
310 + \begin{equation}
311 + T_{jf}=\frac{A_j}{8\pi\eta R_{jf}}\left(I +
312 +  \frac{\mathbf{R}_{jf}\mathbf{R}_{jf}^T}{R_{jf}^2}\right)
313 + \end{equation}
314 + Here, $A_j$ is the area of subfacet $j$ which is a triangle containing
315 + two of the vertices of the facet along with the centroid.
316 + $\mathbf{R}_{jf}$ is the vector between the centroid of facet $f$ and
317 + the centroid of sub-facet $j$, and $I$ is the ($3 \times 3$) identity
318 + matrix.  $\eta$ is the viscosity of the external bath.
319 +
320 + \begin{figure}
321 + \includegraphics[width=\linewidth]{hydro}
322 + \caption{The resistance tensor $\Xi$ for a facet comprising sites $i$,
323 +  $j$, and $k$ is constructed using Oseen tensor contributions between
324 +  the centoid of the facet $f$ and each of the sub-facets ($i,f,j$),
325 +  ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets are
326 +  located at $1$, $2$, and $3$, and the area of each sub-facet is
327 +  easily computed using half the cross product of two of the edges.}
328 + \label{hydro}
329 + \end{figure}
330 +
331 + The Oseen tensors for each of the sub-facets are added together, and
332 + the resulting matrix is inverted to give a $3 \times 3$ resistance
333 + tensor for translations of the triangular facet,
334 + \begin{equation}
335 + \Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}.
336 + \end{equation}
337 + Note that this treatment explicitly ignores rotations (and
338 + translational-rotational coupling) of the facet.  In compact systems,
339 + the facets stay relatively fixed in orientation between
340 + configurations, so this appears to be a reasonably good approximation.
341 +
342 + We have implemented this method by extending the Langevin dynamics
343 + integrator in our code, OpenMD.\cite{Meineke2005,openmd} The Delaunay
344 + triangulation and computation of the convex hull are done using calls
345 + to the qhull library.\cite{qhull} There is a moderate penalty for
346 + computing the convex hull at each step in the molecular dynamics
347 + simulation (HOW MUCH?), but the convex hull is remarkably easy to
348 + parallelize on distributed memory machines (see Appendix A).
349 +
350   \section{Tests \& Applications}
351 + \label{sec:tests}
352  
353 + In order to test this method, we have carried out simulations using
354 + the Langevin Hull on a crystalline system (gold nanoparticles), a
355 + liquid droplet (SPC/E water), and a heterogeneous mixture (gold
356 + nanoparticles in a water droplet).  In each case, we have computed
357 + bulk properties that depend on the external applied pressure.  Of
358 + particular interest is the bulk modulus,
359 + \begin{equation}
360 + \kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right
361 + )_{T}.
362 + \label{eq:BM}
363 + \end{equation}
364 +
365 + One problem with eliminating periodic boundary conditions and
366 + simulation boxes is that the volume of a three-dimensional point cloud
367 + is not well-defined.  In order to compute the compressibility of a
368 + bulk material, we make an assumption that the number density, $\rho =
369 + \frac{N}{V}$, is uniform within some region of the cloud.  The
370 + compressibility can then be expressed in terms of the average number
371 + of particles in that region,
372 + \begin{equation}
373 + \kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right
374 + )_{T}
375 + \label{eq:BMN}
376 + \end{equation}
377 + The region we pick is a spherical volume of 10 \AA radius centered in
378 + the middle of the cluster.  This radius is arbitrary, and any
379 + bulk-like portion of the cluster can be used to compute the bulk
380 + modulus.
381 +
382 + One might also assume that the volume of the convex hull could be
383 + taken as the system volume in the compressibility expression
384 + (Eq. \ref{eq:BM}), but this has implications at lower pressures (which
385 + are explored in detail in the section on water droplets).
386 +
387   \subsection{Bulk modulus of gold nanoparticles}
388  
389   \begin{figure}
# Line 199 | Line 413 | Both NVT \cite{Glattli2002} and NPT \cite{Motakabbir19
413  
414   \subsection{Compressibility of SPC/E water clusters}
415  
416 < 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.
416 > Prior molecular dynamics simulations on SPC/E water (both in
417 > NVT~\cite{Glattli2002} and NPT~\cite{Motakabbir1990, Pi2009}
418 > ensembles) have yielded values for the isothermal compressibility that
419 > agree well with experiment.\cite{Fine1973} The results of two
420 > different approaches for computing the isothermal compressibility from
421 > Langevin Hull simulations for pressures between 1 and 6500 atm are
422 > shown in Fig. \ref{fig:compWater} along with compressibility values
423 > obtained from both other SPC/E simulations and experiment.
424 > Compressibility values from all references are for applied pressures
425 > within the range 1 - 1000 atm.
426  
427   \begin{figure}
428 < \includegraphics[width=\linewidth]{new_isothermal}
428 > \includegraphics[width=\linewidth]{new_isothermalN}
429   \caption{Compressibility of SPC/E water}
430 < \label{compWater}
430 > \label{fig:compWater}
431   \end{figure}
432  
433 < We initially used the classic compressibility formula
433 > Isothermal compressibility values calculated using the number density
434 > (Eq. \ref{eq:BMN}) expression are in good agreement with experimental
435 > and previous simulation work throughout the 1 - 1000 atm pressure
436 > regime.  Compressibilities computed using the Hull volume, however,
437 > deviate dramatically from the experimental values at low applied
438 > pressures.  The reason for this deviation is quite simple; at low
439 > applied pressures, the liquid is in equilibrium with a vapor phase,
440 > and it is entirely possible for one (or a few) molecules to drift away
441 > from the liquid cluster (see Fig. \ref{fig:coneOfShame}).  At low
442 > pressures, the restoring forces on the facets are very gentle, and
443 > this means that the hulls often take on relatively distorted
444 > geometries which include large volumes of empty space.
445  
446 < \begin{equation}
447 < \kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T}
448 < \end{equation}
446 > \begin{figure}
447 > \includegraphics[width=\linewidth]{flytest2}
448 > \caption{At low pressures, the liquid is in equilibrium with the vapor
449 >  phase, and isolated molecules can detach from the liquid droplet.
450 >  This is expected behavior, but the reported volume of the convex
451 >  hull includes large regions of empty space.  For this reason,
452 >  compressibilities should be computed using local number densities
453 >  rather than hull volumes.}
454 > \label{fig:coneOfShame}
455 > \end{figure}
456  
457 < 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.
457 > At higher pressures, the equilibrium favors the liquid phase, and the
458 > hull geometries are much more compact.  Because of the liquid-vapor
459 > effect on the convex hull, the regional number density approach
460 > (Eq. \ref{eq:BMN}) provides more reliable estimates of the bulk
461 > modulus.
462  
463 < Per the fluctuation dissipation theorem \cite{Debendedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure
463 > 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.
464  
465 + Regardless of the difficulty in obtaining accurate hull
466 + volumes at low temperature and pressures, the Langevin Hull NPT method
467 + provides reasonable isothermal compressibility values for water
468 + through a large range of pressures.
469 +
470 + 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
471 +
472   \begin{equation}
473   \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T}
474   \end{equation}
475  
476   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.
477  
226 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
478  
228 \begin{equation}
229 \kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T}
230 \end{equation}
231
232 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.
233
479   \subsection{Molecular orientation distribution at cluster boundary}
480  
481   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.
# Line 275 | Line 520 | The orientational preference exhibited by hull molecul
520  
521   \section{Appendix A: Hydrodynamic tensor for triangular facets}
522  
278 \begin{figure}
279 \includegraphics[width=\linewidth]{hydro}
280 \caption{Hydro}
281 \label{hydro}
282 \end{figure}
283
284 \begin{equation}
285 \Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}
286 \end{equation}
287
288 \begin{equation}
289 T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I +
290  \frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right)
291 \end{equation}
292
523   \section{Appendix B: Computing Convex Hulls on Parallel Computers}
524  
525   \section{Acknowledgments}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines