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 3665 by gezelter, Wed Oct 27 18:48:34 2010 UTC vs.
Revision 3689 by gezelter, Thu Nov 11 22:17:37 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 40 | Line 40 | Notre Dame, Indiana 46556}
40    We have developed a new isobaric-isothermal (NPT) algorithm which
41    applies an external pressure to the facets comprising the convex
42    hull surrounding the system.  A Langevin thermostat is also applied
43 <  to facets of the hull to mimic contact with an external heat
44 <  bath. This new method, the ``Langevin Hull'', performs better than
45 <  traditional affine transform methods for systems containing
46 <  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 <  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.
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 maintain a ta
66   \section{Introduction}
67  
68   The most common molecular dynamics methods for sampling configurations
69 < of an isobaric-isothermal (NPT) ensemble maintain a target pressure in
70 < a simulation by coupling the volume of the system to a {\it barostat},
71 < which is an extra degree of freedom propagated along with the particle
72 < coordinates.  These methods require periodic boundary conditions,
73 < because when the instantaneous pressure in the system differs from the
74 < target pressure, the volume is reduced or expanded using {\it affine
75 <  transforms} of the system geometry. An affine transform scales the
76 < size and shape of the periodic box as well as the particle positions
77 < within the box (but not the sizes of the particles). The most common
78 < constant pressure methods, including the Melchionna
79 < modification\cite{Melchionna1993} to the Nos\'e-Hoover-Andersen
80 < equations of motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx}
81 < the Berendsen pressure bath,\cite{ISI:A1984TQ73500045} and the
82 < Langevin Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize
69 > from an isobaric-isothermal (NPT) ensemble maintain a target pressure
70 > in a simulation by coupling the volume of the system to a {\it
71 >  barostat}, which is an extra degree of freedom propagated along with
72 > the particle coordinates.  These methods require periodic boundary
73 > conditions, because when the instantaneous pressure in the system
74 > differs from the target pressure, the volume is reduced or expanded
75 > using {\it affine transforms} of the system geometry. An affine
76 > transform scales the size and shape of the periodic box as well as the
77 > particle positions within the box (but not the sizes of the
78 > particles). The most common constant pressure methods, including the
79 > Melchionna modification\cite{Melchionna1993} to the
80 > Nos\'e-Hoover-Andersen equations of
81 > motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen
82 > pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin
83 > Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize scaled
84   coordinate transformation to adjust the box volume.  As long as the
85 < material in the simulation box is essentially a bulk-like liquid which
86 < has a relatively uniform compressibility, the standard affine
87 < transform approach provides an excellent way of adjusting the volume
88 < of the system and applying pressure directly via the interactions
88 < between atomic sites.
85 > material in the simulation box has a relatively uniform
86 > compressibility, the standard affine transform approach provides an
87 > excellent way of adjusting the volume of the system and applying
88 > pressure directly via the interactions between atomic sites.
89  
90   One problem with this approach appears when the system being simulated
91   is an inhomogeneous mixture in which portions of the simulation box
# Line 100 | Line 100 | slow enough to avoid the instabilities in the incompre
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  
# Line 115 | Line 115 | effect.  For example, calculations using typical hydra
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 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 + Typically {\it total} protein concentrations in the cell are on the
125 + order of 160-310 mg/ml,\cite{Brown1991195} and individual proteins
126 + have concentrations orders of magnitude lower than this in the
127 + cellular environment. The effective concentrations of single proteins
128 + in simulations may have significant effects on the structure and
129 + dynamics of simulated structures.
130 +
131   \subsection*{Boundary Methods}
132 < There have been a number of other approaches to explicit
133 < non-periodicity that focus on constant or nearly-constant {\it volume}
134 < conditions while maintaining bulk-like behavior.  Berkowitz and
135 < McCammon introduced a stochastic (Langevin) boundary layer inside a
136 < region of fixed molecules which effectively enforces constant
137 < temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this
138 < approach, the stochastic and fixed regions were defined relative to a
139 < central atom.  Brooks and Karplus extended this method to include
140 < deformable stochastic boundaries.\cite{iii:6312} The stochastic
141 < boundary approach has been used widely for protein
142 < simulations. [CITATIONS NEEDED]
132 > There have been a number of approaches to handle simulations of
133 > explicitly non-periodic systems that focus on constant or
134 > nearly-constant {\it volume} conditions while maintaining bulk-like
135 > behavior.  Berkowitz and McCammon introduced a stochastic (Langevin)
136 > boundary layer inside a region of fixed molecules which effectively
137 > enforces constant temperature and volume (NVT)
138 > conditions.\cite{Berkowitz1982} In this approach, the stochastic and
139 > fixed regions were defined relative to a central atom.  Brooks and
140 > Karplus extended this method to include deformable stochastic
141 > boundaries.\cite{iii:6312} The stochastic boundary approach has been
142 > used widely for protein simulations.
143  
144   The electrostatic and dispersive behavior near the boundary has long
145   been a cause for concern when performing simulations of explicitly
# Line 142 | Line 150 | affecting most of molecules in the simulation.  This r
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 molecules in the simulation.  This result echoes an
154 < earlier study which showed that an extended planar hydrophobic surface
155 < caused orientational preference at the surface which extended
156 < relatively deep (7 \r{A}) into the liquid simulation
157 < cell.\cite{Lee1984} The surface constrained all-atom solvent (SCAAS)
158 < model \cite{King1989} improved upon its SCSSD predecessor. The SCAAS
159 < model utilizes a polarization constraint which is applied to the
160 < surface molecules to maintain bulk-like structure at the cluster
161 < surface. A radial constraint is used to maintain the desired bulk
162 < density of the liquid. Both constraint forces are applied only to a
163 < pre-determined number of the outermost molecules.
153 > affecting most of the molecules in the simulation.  This result echoes
154 > an earlier study which showed that an extended planar hydrophobic
155 > surface caused orientational preferences at the surface which extended
156 > relatively deep (7 \AA) into the liquid simulation cell.\cite{Lee1984}
157 > The surface constrained all-atom solvent (SCAAS) model \cite{King1989}
158 > improved upon its SCSSD predecessor. The SCAAS model utilizes a
159 > polarization constraint which is applied to the surface molecules to
160 > maintain bulk-like structure at the cluster surface. A radial
161 > constraint is used to maintain the desired bulk density of the
162 > liquid. Both constraint forces are applied only to a pre-determined
163 > number of the outermost molecules.
164  
165   Beglov and Roux have developed a boundary model in which the hard
166   sphere boundary has a radius that varies with the instantaneous
167   configuration of the solute (and solvent) molecules.\cite{beglov:9050}
168   This model contains a clear pressure and surface tension contribution
169 < to the free energy which XXX.
169 > to the free energy.
170  
171   \subsection*{Restraining Potentials}
172   Restraining {\it potentials} introduce repulsive potentials at the
# Line 173 | Line 181 | position of the nearest solute atom.\cite{LiY._jp04685
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:xw} This
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
# Line 189 | Line 197 | force in a direction that is inward-facing relative to
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 202 | 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.  The methodology is introduced in section \ref{sec:meth}, tests
217 < on crystalline nanoparticles, liquid clusters, and heterogeneous
218 < mixtures are detailed in section \ref{sec:tests}.  Section
219 < \ref{sec:discussion} summarizes our findings.
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$).  This bath interacts only with the
226 < objects on the exterior hull of the system.  Defining the hull of the
227 < simulation is done in a manner similar to the approach of Kohanoff,
228 < Caro and Finnis.\cite{Kohanoff:2005qm} That is, any instantaneous
229 < configuration of the atoms in the system is considered as a point
230 < cloud in three dimensional space.  Delaunay triangulation is used to
231 < find all facets between coplanar
232 < neighbors.\cite{delaunay,springerlink:10.1007/BF00977785}  In highly
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 3-space
236 < that contain exactly three atoms.
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
# Line 238 | Line 247 | simulation.
247   simulation.
248  
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  
# Line 293 | Line 302 | viscosity of the fluid.  The resistance tensor is rela
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 323 | 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 363 | Line 372 | Note that this treatment explicitly ignores rotations
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.
# Line 373 | Line 382 | molecular dynamics time step, the following process is
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 done using the current atomic
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:
# Line 447 | Line 456 | atoms and the SPC/E water molecules.\cite{ISI:00016776
456   Spohr potential was adopted in depicting the interaction between metal
457   atoms and the SPC/E water molecules.\cite{ISI:000167766600035}
458  
459 < \subsection{Bulk modulus of gold nanoparticles}
459 > \subsection{Bulk Modulus of gold nanoparticles}
460  
461 < The compressibility is well-known for gold, and it provides a good first
462 < test of how the method compares to other similar methods.  
463 <
464 < \begin{figure}
465 < \includegraphics[width=\linewidth]{P_T_combined}
466 < \caption{Pressure and temperature response of an 18 \AA\ gold
467 <  nanoparticle initially when first placed in the Langevin Hull
459 <  ($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa) and starting
460 <  from initial conditions that were far from the bath pressure and
461 <  temperature.  The pressure response is rapid, and the thermal
462 <  equilibration depends on both total surface area and the viscosity
463 <  of the bath.}
464 < \label{pressureResponse}
465 < \end{figure}
461 > The compressibility (and its inverse, the bulk modulus) is well-known
462 > for gold, and is captured well by the embedded atom method
463 > (EAM)~\cite{PhysRevB.33.7983} potential and related multi-body force
464 > fields.  In particular, the quantum Sutton-Chen potential gets nearly
465 > quantitative agreement with the experimental bulk modulus values, and
466 > makes a good first test of how the Langevin Hull will perform at large
467 > applied pressures.
468  
469 + The Sutton-Chen (SC) potentials are based on a model of a metal which
470 + treats the nuclei and core electrons as pseudo-atoms embedded in the
471 + electron density due to the valence electrons on all of the other
472 + atoms in the system.\cite{Chen90} The SC potential has a simple form
473 + that closely resembles the Lennard Jones potential,
474   \begin{equation}
475 < \kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial
476 <    P}\right)
475 > \label{eq:SCP1}
476 > 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] ,
477   \end{equation}
478 + where $V^{pair}_{ij}$ and $\rho_{i}$ are given by
479 + \begin{equation}
480 + \label{eq:SCP2}
481 + 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}}.
482 + \end{equation}
483 + $V^{pair}_{ij}$ is a repulsive pairwise potential that accounts for
484 + interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in
485 + Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models
486 + the interactions between the valence electrons and the cores of the
487 + pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy
488 + scale, $c_i$ scales the attractive portion of the potential relative
489 + to the repulsive interaction and $\alpha_{ij}$ is a length parameter
490 + that assures a dimensionless form for $\rho$. These parameters are
491 + tuned to various experimental properties such as the density, cohesive
492 + energy, and elastic moduli for FCC transition metals. The quantum
493 + Sutton-Chen (QSC) formulation matches these properties while including
494 + zero-point quantum corrections for different transition
495 + metals.\cite{PhysRevB.59.3527,QSC}
496  
497 + In bulk gold, the experimentally-measured value for the bulk modulus
498 + is 180.32 GPa, while previous calculations on the QSC potential in
499 + periodic-boundary simulations of the bulk crystal have yielded values
500 + of 175.53 GPa.\cite{QSC} Using the same force field, we have performed
501 + a series of relatively short (200 ps) simulations on 40 \AA~ radius
502 + nanoparticles under the Langevin Hull at a variety of applied
503 + pressures ranging from 0 -- 10 GPa.  We obtain a value of 177.55 GPa
504 + for the bulk modulus of gold using this technique, in close agreement
505 + with both previous simulations and the experimental bulk modulus of
506 + gold.
507 +
508   \begin{figure}
509 < \includegraphics[width=\linewidth]{compress_tb}
510 < \caption{Isothermal Compressibility (18 \AA gold nanoparticle)}
511 < \label{temperatureResponse}
509 > \includegraphics[width=\linewidth]{stacked}
510 > \caption{The response of the internal pressure and temperature of gold
511 >  nanoparticles when first placed in the Langevin Hull
512 >  ($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa), starting
513 >  from initial conditions that were far from the bath pressure and
514 >  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).}
515 > \label{fig:pressureResponse}
516   \end{figure}
517 +
518 + We note that the Langevin Hull produces rapidly-converging behavior
519 + for structures that are started far from equilibrium.  In
520 + Fig. \ref{fig:pressureResponse} we show how the pressure and
521 + temperature respond to the Langevin Hull for nanoparticles that were
522 + initialized far from the target pressure and temperature.  As
523 + expected, the rate at which thermal equilibrium is achieved depends on
524 + the total surface area of the cluter exposed to the bath as well as
525 + the bath viscosity.  Pressure that is applied suddenly to a cluster
526 + can excite breathing vibrations, but these rapidly damp out (on time
527 + scales of 30-50 ps).
528  
529   \subsection{Compressibility of SPC/E water clusters}
530  
# Line 485 | Line 536 | Compressibility values from all references are for app
536   Langevin Hull simulations for pressures between 1 and 6500 atm are
537   shown in Fig. \ref{fig:compWater} along with compressibility values
538   obtained from both other SPC/E simulations and experiment.
488 Compressibility values from all references are for applied pressures
489 within the range 1 - 1000 atm.
539  
540   \begin{figure}
541   \includegraphics[width=\linewidth]{new_isothermalN}
# Line 496 | Line 545 | and previous simulation work throughout the 1 - 1000 a
545  
546   Isothermal compressibility values calculated using the number density
547   (Eq. \ref{eq:BMN}) expression are in good agreement with experimental
548 < and previous simulation work throughout the 1 - 1000 atm pressure
548 > and previous simulation work throughout the 1 -- 1000 atm pressure
549   regime.  Compressibilities computed using the Hull volume, however,
550   deviate dramatically from the experimental values at low applied
551   pressures.  The reason for this deviation is quite simple; at low
# Line 508 | Line 557 | geometries which include large volumes of empty space.
557   geometries which include large volumes of empty space.
558  
559   \begin{figure}
560 < \includegraphics[width=\linewidth]{flytest2}
560 > \includegraphics[width=\linewidth]{coneOfShame}
561   \caption{At low pressures, the liquid is in equilibrium with the vapor
562    phase, and isolated molecules can detach from the liquid droplet.
563    This is expected behavior, but the volume of the convex hull
564 <  includes large regions of empty space.  For this reason,
564 >  includes large regions of empty space. For this reason,
565    compressibilities are computed using local number densities rather
566    than hull volumes.}
567   \label{fig:coneOfShame}
# Line 522 | Line 571 | bulk modulus.
571   and the hull geometries are much more compact.  Because of the
572   liquid-vapor effect on the convex hull, the regional number density
573   approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the
574 < bulk modulus.
574 > compressibility.
575  
576   In both the traditional compressibility formula (Eq. \ref{eq:BM}) and
577   the number density version (Eq. \ref{eq:BMN}), multiple simulations at
# Line 533 | Line 582 | volume,\cite{Debenedetti1986},
582   \begin{equation}
583   \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle
584      V \right \rangle ^{2}}{V \, k_{B} \, T},
585 + \label{eq:BMVfluct}
586   \end{equation}
587   or, equivalently, fluctuations in the number of molecules within the
588   fixed region,
589   \begin{equation}
590   \kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle
591      N \right \rangle ^{2}}{N \, k_{B} \, T},
592 + \label{eq:BMNfluct}
593   \end{equation}
594   Thus, the compressibility of each simulation can be calculated
595 < entirely independently from all other trajectories. However, the
596 < resulting compressibilities were still as much as an order of
597 < magnitude larger than the reference values.  Any compressibility
598 < calculation that relies on the hull volume will suffer these effects.
599 < WE NEED MORE HERE.
595 > entirely independently from other trajectories.  Compressibility
596 > calculations that rely on the hull volume will still suffer the
597 > effects of the empty space due to the vapor phase; for this reason, we
598 > recommend using the number density (Eq. \ref{eq:BMN}) or number
599 > density fluctuations (Eq. \ref{eq:BMNfluct}) for computing
600 > compressibilities.
601  
602   \subsection{Molecular orientation distribution at cluster boundary}
603  
604 < In order for non-periodic boundary conditions to be widely applicable,
604 > In order for a non-periodic boundary method to be widely applicable,
605   they must be constructed in such a way that they allow a finite system
606 < to replicate the properties of the bulk.  Naturally, this requirement
607 < has spawned many methods for fixing and characterizing the effects of
608 < artifical boundaries. Of particular interest regarding the Langevin
609 < Hull is the orientation of water molecules that are part of the
610 < geometric hull.  Ideally, all molecules in the cluster will have the
611 < same orientational distribution as bulk water.
606 > to replicate the properties of the bulk. Early non-periodic simulation
607 > methods (e.g. hydrophobic boundary potentials) induced spurious
608 > orientational correlations deep within the simulated
609 > system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
610 > fixing and characterizing the effects of artifical boundaries
611 > including methods which fix the orientations of a set of edge
612 > molecules.\cite{Warshel1978,King1989}
613  
614 < The orientation of molecules at the edges of a simulated cluster has
615 < long been a concern when performing simulations of explicitly
616 < non-periodic systems. Early work led to the surface constrained soft
617 < sphere dipole model (SCSSD) \cite{Warshel1978} in which the surface
618 < molecules are fixed in a random orientation representative of the bulk
619 < solvent structural properties. Belch, et al \cite{Belch1985} simulated
620 < clusters of TIPS2 water surrounded by a hydrophobic bounding
621 < potential. The spherical hydrophobic boundary induced dangling
622 < hydrogen bonds at the surface that propagated deep into the cluster,
623 < affecting 70\% of the 100 molecules in the simulation. This result
624 < echoes an earlier study which showed that an extended planar
572 < hydrophobic surface caused orientational preference at the surface
573 < which extended 7 \r{A} into the liquid simulation cell
574 < \cite{Lee1984}. The surface constrained all-atom solvent (SCAAS) model
575 < \cite{King1989} improved upon its SCSSD predecessor. The SCAAS model
576 < utilizes a polarization constraint which is applied to the surface
577 < molecules to maintain bulk-like structure at the cluster surface. A
578 < radial constraint is used to maintain the desired bulk density of the
579 < liquid. Both constraint forces are applied only to a pre-determined
580 < number of the outermost molecules.
614 > As described above, the Langevin Hull does not require that the
615 > orientation of molecules be fixed, nor does it utilize an explicitly
616 > hydrophobic boundary, or orientational or radial constraints.
617 > Therefore, the orientational correlations of the molecules in water
618 > clusters are of particular interest in testing this method.  Ideally,
619 > the water molecules on the surfaces of the clusterss will have enough
620 > mobility into and out of the center of the cluster to maintain
621 > bulk-like orientational distribution in the absence of orientational
622 > and radial constraints.  However, since the number of hydrogen bonding
623 > partners available to molecules on the exterior are limited, it is
624 > likely that there will be an effective hydrophobicity of the hull.
625  
626 < In contrast, the Langevin Hull does not require that the orientation
627 < of molecules be fixed, nor does it utilize an explicitly hydrophobic
628 < boundary, orientational constraint or radial constraint. The number
629 < and identity of the molecules included on the convex hull are dynamic
586 < properties, thus avoiding the formation of an artificial solvent
587 < boundary layer. The hope is that the water molecules on the surface of
588 < the cluster, if left to their own devices in the absence of
589 < orientational and radial constraints, will maintain a bulk-like
590 < orientational distribution.
591 <
592 < 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.
593 <
594 < The orientation of a water molecule is described by
595 <
626 > To determine the extent of these effects, we examined the
627 > orientationations exhibited by SPC/E water in a cluster of 1372
628 > molecules at 300 K and at pressures ranging from 1 -- 1000 atm.  The
629 > orientational angle of a water molecule is described
630   \begin{equation}
631   \cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
632   \end{equation}
633 + where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of
634 + mass and the cluster center of mass, and $\vec{\mu}_{i}$ is the vector
635 + bisecting the H-O-H angle of molecule {\it i}.  Bulk-like
636 + distributions will result in $\langle \cos \theta \rangle$ values
637 + close to zero.  If the hull exhibits an overabundance of
638 + externally-oriented oxygen sites, the average orientation will be
639 + negative, while dangling hydrogen sites will result in positive
640 + average orientations.
641  
642 < 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}.
643 <
642 > Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
643 > for molecules in the interior of the cluster (squares) and for
644 > molecules included in the convex hull (circles).
645   \begin{figure}
603 \includegraphics[width=\linewidth]{g_r_theta}
604 \caption{Definition of coordinates}
605 \label{coords}
606 \end{figure}
607
608 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).
609
610 \begin{figure}
646   \includegraphics[width=\linewidth]{pAngle}
647 < \caption{SPC/E water clusters: only minor dewetting at the boundary}
648 < \label{pAngle}
647 > \caption{Distribution of $\cos{\theta}$ values for molecules on the
648 >  interior of the cluster (squares) and for those participating in the
649 >  convex hull (circles) at a variety of pressures.  The Langevin hull
650 >  exhibits minor dewetting behavior with exposed oxygen sites on the
651 >  hull water molecules.  The orientational preference for exposed
652 >  oxygen appears to be independent of applied pressure. }
653 > \label{fig:pAngle}
654   \end{figure}
655  
656 < 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.
656 > As expected, interior molecules (those not included in the convex
657 > hull) maintain a bulk-like structure with a uniform distribution of
658 > orientations. Molecules included in the convex hull show a slight
659 > preference for values of $\cos{\theta} < 0.$ These values correspond
660 > to molecules with oxygen directed toward the exterior of the cluster,
661 > forming a dangling hydrogen bond acceptor site.
662  
663 < 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.
663 > In the absence of an electrostatic contribution from the exterior
664 > bath, the orientational distribution of water molecules included in
665 > the Langevin Hull will slightly resemble the distribution at a neat
666 > water liquid/vapor interface.  Previous molecular dynamics simulations
667 > of SPC/E water \cite{Taylor1996} have shown that molecules at the
668 > liquid/vapor interface favor an orientation where one hydrogen
669 > protrudes from the liquid phase. This behavior is demonstrated by
670 > experiments \cite{Du1994} \cite{Scatena2001} showing that
671 > approximately one-quarter of water molecules at the liquid/vapor
672 > interface form dangling hydrogen bonds. The negligible preference
673 > shown in these cluster simulations could be removed through the
674 > introduction of an implicit solvent model, which would provide the
675 > missing electrostatic interactions between the cluster molecules and
676 > the surrounding temperature/pressure bath.
677  
678 < 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.
678 > The orientational preference exhibited by hull molecules in the
679 > Langevin hull is significantly weaker than the preference caused by an
680 > explicit hydrophobic bounding potential.  Additionally, the Langevin
681 > Hull does not require that the orientation of any molecules be fixed
682 > in order to maintain bulk-like structure, even at the cluster surface.
683  
684   \subsection{Heterogeneous nanoparticle / water mixtures}
685  
686 + To further test the method, we simulated gold nanopartices ($r = 18$
687 + \AA) solvated by explicit SPC/E water clusters using the Langevin
688 + hull.  This was done at pressures of 1, 2, 5, 10, 20, 50 and 100 atm
689 + in order to observe the effects of pressure on the ordering of water
690 + ordering at the surface.  In Fig. \ref{fig:RhoR} we show the density
691 + of water adjacent to the surface as a function of pressure, as well as
692 + the orientational ordering of water at the surface of the
693 + nanoparticle.
694 +
695 + \begin{figure}
696 +
697 + \caption{interesting plot showing cluster behavior}
698 + \label{fig:RhoR}
699 + \end{figure}
700 +
701 + At higher pressures, problems with the gold - water interaction
702 + potential became apparent.  The model we are using (due to Spohr) was
703 + intended for relatively low pressures; it utilizes both shifted Morse
704 + and repulsive Morse potentials to model the Au/O and Au/H
705 + interactions, respectively.  The repulsive wall of the Morse potential
706 + does not diverge quickly enough at short distances to prevent water
707 + from diffusing into the center of the gold nanoparticles.  This
708 + behavior is likely not a realistic description of the real physics of
709 + the situation.  A better model of the gold-water adsorption behavior
710 + appears to require harder repulsive walls to prevent this behavior.
711 +
712   \section{Discussion}
713   \label{sec:discussion}
714  
715 + The Langevin Hull samples the isobaric-isothermal ensemble for
716 + non-periodic systems by coupling the system to a bath characterized by
717 + pressure, temperature, and solvent viscosity.  This enables the
718 + simulation of heterogeneous systems composed of materials of
719 + significantly different compressibilities.  Because the boundary is
720 + dynamically determined during the simulation and the molecules
721 + interacting with the boundary can change, the method and has minimal
722 + perturbations on the behavior of molecules at the edges of the
723 + simulation.  Further work on this method will involve implicit
724 + electrostatics at the boundary (which is missing in the current
725 + implementation) as well as more sophisticated treatments of the
726 + surface geometry (alpha
727 + shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight
728 + Cocone\cite{Dey:2003ts}). The non-convex hull geometries are
729 + significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull
730 + ($\mathcal{O}(N \log N)$), but would enable the use of hull volumes
731 + directly in computing the compressibility of the sample.
732 +
733   \section*{Appendix A: Computing Convex Hulls on Parallel Computers}
734  
735 + In order to use the Langevin Hull for simulations on parallel
736 + computers, one of the more difficult tasks is to compute the bounding
737 + surface, facets, and resistance tensors when the individual processors
738 + have incomplete information about the entire system's topology.  Most
739 + parallel decomposition methods assign primary responsibility for the
740 + motion of an atomic site to a single processor, and we can exploit
741 + this to efficiently compute the convex hull for the entire system.
742 +
743 + The basic idea involves splitting the point cloud into
744 + spatially-overlapping subsets and computing the convex hulls for each
745 + of the subsets.  The points on the convex hull of the entire system
746 + are all present on at least one of the subset hulls. The algorithm
747 + works as follows:
748 + \begin{enumerate}
749 + \item Each processor computes the convex hull for its own atomic sites
750 +  (left panel in Fig. \ref{fig:parallel}).
751 + \item The Hull vertices from each processor are communicated to all of
752 +  the processors, and each processor assembles a complete list of hull
753 +  sites (this is much smaller than the original number of points in
754 +  the point cloud).
755 + \item Each processor computes the global convex hull (right panel in
756 +  Fig. \ref{fig:parallel}) using only those points that are the union
757 +  of sites gathered from all of the subset hulls.  Delaunay
758 +  triangulation is then done to obtain the facets of the global hull.
759 + \end{enumerate}
760 +
761 + \begin{figure}
762 + \includegraphics[width=\linewidth]{parallel}
763 + \caption{When the sites are distributed among many nodes for parallel
764 +  computation, the processors first compute the convex hulls for their
765 +  own sites (dashed lines in left panel). The positions of the sites
766 +  that make up the subset hulls are then communicated to all
767 +  processors (middle panel).  The convex hull of the system (solid line in
768 +  right panel) is the convex hull of the points on the union of the subset
769 +  hulls.}
770 + \label{fig:parallel}
771 + \end{figure}
772 +
773 + The individual hull operations scale with
774 + $\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total
775 + number of sites, and $p$ is the number of processors.  These local
776 + hull operations create a set of $p$ hulls each with approximately
777 + $\frac{n}{3pr}$ sites (for a cluster of radius $r$). The worst-case
778 + communication cost for using a ``gather'' operation to distribute this
779 + information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n
780 +  \beta (p-1)}{3 r p^2})$, while the final computation of the system
781 + hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$.
782 +
783 + For a large number of atoms on a moderately parallel machine, the
784 + total costs are dominated by the computations of the individual hulls,
785 + and communication of these hulls to so the Langevin hull sees roughly
786 + linear speed-up with increasing processor counts.
787 +
788   \section*{Acknowledgments}
789   Support for this project was provided by the
790   National Science Foundation under grant CHE-0848243. Computational
791   time was provided by the Center for Research Computing (CRC) at the
792   University of Notre Dame.  
793  
794 + Molecular graphics images were produced using the UCSF Chimera package from
795 + the Resource for Biocomputing, Visualization, and Informatics at the
796 + University of California, San Francisco (supported by NIH P41 RR001081).
797   \newpage
798  
799   \bibliography{langevinHull}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines