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 3668 by gezelter, Thu Nov 4 16:01:58 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 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 121 | Line 121 | There have been a number of other approaches to explic
121   SIMULATIONS]
122  
123   \subsection*{Boundary Methods}
124 < There have been a number of other approaches to explicit
125 < non-periodicity that focus on constant or nearly-constant {\it volume}
126 < conditions while maintaining bulk-like behavior.  Berkowitz and
127 < McCammon introduced a stochastic (Langevin) boundary layer inside a
128 < region of fixed molecules which effectively enforces constant
129 < temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this
130 < approach, the stochastic and fixed regions were defined relative to a
131 < central atom.  Brooks and Karplus extended this method to include
132 < deformable stochastic boundaries.\cite{iii:6312} The stochastic
133 < boundary approach has been used widely for protein
134 < simulations. [CITATIONS NEEDED]
124 > There have been a number of approaches to handle simulations of
125 > explicitly non-periodic systems that focus on constant or
126 > nearly-constant {\it volume} conditions while maintaining bulk-like
127 > behavior.  Berkowitz and McCammon introduced a stochastic (Langevin)
128 > boundary layer inside a region of fixed molecules which effectively
129 > enforces constant temperature and volume (NVT)
130 > conditions.\cite{Berkowitz1982} In this approach, the stochastic and
131 > fixed regions were defined relative to a central atom.  Brooks and
132 > Karplus extended this method to include deformable stochastic
133 > boundaries.\cite{iii:6312} The stochastic boundary approach has been
134 > used widely for protein simulations. [CITATIONS NEEDED]
135  
136   The electrostatic and dispersive behavior near the boundary has long
137   been a cause for concern when performing simulations of explicitly
# Line 202 | Line 202 | random forces on the facets of the {\it hull itself} i
202   In the following sections, we extend and generalize the approach of
203   Kohanoff, Caro, and Finnis. The new method, which we are calling the
204   ``Langevin Hull'' applies the external pressure, Langevin drag, and
205 < random forces on the facets of the {\it hull itself} instead of the
206 < atomic sites comprising the vertices of the hull.  This allows us to
207 < decouple the external pressure contribution from the drag and random
208 < force.  The methodology is introduced in section \ref{sec:meth}, tests
209 < on crystalline nanoparticles, liquid clusters, and heterogeneous
210 < mixtures are detailed in section \ref{sec:tests}.  Section
211 < \ref{sec:discussion} summarizes our findings.
205 > random forces on the {\it facets of the hull} instead of the atomic
206 > sites comprising the vertices of the hull.  This allows us to decouple
207 > the external pressure contribution from the drag and random force.
208 > The methodology is introduced in section \ref{sec:meth}, tests on
209 > crystalline nanoparticles, liquid clusters, and heterogeneous mixtures
210 > are detailed in section \ref{sec:tests}.  Section \ref{sec:discussion}
211 > summarizes our findings.
212  
213   \section{Methodology}
214   \label{sec:meth}
# Line 241 | Line 241 | simulation.
241   \includegraphics[width=\linewidth]{hullSample}
242   \caption{The external temperature and pressure bath interacts only
243    with those atoms on the convex hull (grey surface).  The hull is
244 <  computed dynamically at each time step, and molecules dynamically
245 <  move between the interior (Newtonian) region and the Langevin hull.}
244 >  computed dynamically at each time step, and molecules can move
245 >  between the interior (Newtonian) region and the Langevin hull.}
246   \label{fig:hullSample}
247   \end{figure}
248  
# Line 363 | Line 363 | Note that this treatment explicitly ignores rotations
363   \begin{equation}
364   \Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}.
365   \end{equation}
366 < Note that this treatment explicitly ignores rotations (and
366 > Note that this treatment ignores rotations (and
367   translational-rotational coupling) of the facet.  In compact systems,
368   the facets stay relatively fixed in orientation between
369   configurations, so this appears to be a reasonably good approximation.
# Line 373 | Line 373 | molecular dynamics time step, the following process is
373   molecular dynamics time step, the following process is carried out:
374   \begin{enumerate}
375   \item The standard inter-atomic forces ($\nabla_iU$) are computed.
376 < \item Delaunay triangulation is done using the current atomic
376 > \item Delaunay triangulation is carried out using the current atomic
377    configuration.
378   \item The convex hull is computed and facets are identified.
379   \item For each facet:
# Line 447 | Line 447 | atoms and the SPC/E water molecules.\cite{ISI:00016776
447   Spohr potential was adopted in depicting the interaction between metal
448   atoms and the SPC/E water molecules.\cite{ISI:000167766600035}
449  
450 < \subsection{Bulk modulus of gold nanoparticles}
450 > \subsection{Compressibility of gold nanoparticles}
451  
452   The compressibility is well-known for gold, and it provides a good first
453   test of how the method compares to other similar methods.  
# Line 522 | Line 522 | bulk modulus.
522   and the hull geometries are much more compact.  Because of the
523   liquid-vapor effect on the convex hull, the regional number density
524   approach (Eq. \ref{eq:BMN}) provides more reliable estimates of the
525 < bulk modulus.
525 > compressibility.
526  
527   In both the traditional compressibility formula (Eq. \ref{eq:BM}) and
528   the number density version (Eq. \ref{eq:BMN}), multiple simulations at
# Line 543 | Line 543 | magnitude larger than the reference values.  Any compr
543   Thus, the compressibility of each simulation can be calculated
544   entirely independently from all other trajectories. However, the
545   resulting compressibilities were still as much as an order of
546 < magnitude larger than the reference values.  Any compressibility
546 > magnitude larger than the reference values. However, compressibility
547   calculation that relies on the hull volume will suffer these effects.
548   WE NEED MORE HERE.
549  
# Line 551 | Line 551 | to replicate the properties of the bulk.  Naturally, t
551  
552   In order for non-periodic boundary conditions to be widely applicable,
553   they must be constructed in such a way that they allow a finite system
554 < to replicate the properties of the bulk.  Naturally, this requirement
555 < has spawned many methods for fixing and characterizing the effects of
556 < artifical boundaries. Of particular interest regarding the Langevin
557 < Hull is the orientation of water molecules that are part of the
558 < geometric hull.  Ideally, all molecules in the cluster will have the
559 < same orientational distribution as bulk water.
554 > to replicate the properties of the bulk.  Early non-periodic
555 > simulation methods (e.g. hydrophobic boundary potentials) induced
556 > spurious orientational correlations deep within the simulated
557 > system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
558 > fixing and characterizing the effects of artifical boundaries
559 > including methods which fix the orientations of a set of edge
560 > molecules.\cite{Warshel1978,King1989}
561  
562 < The orientation of molecules at the edges of a simulated cluster has
563 < long been a concern when performing simulations of explicitly
564 < non-periodic systems. Early work led to the surface constrained soft
565 < sphere dipole model (SCSSD) \cite{Warshel1978} in which the surface
566 < molecules are fixed in a random orientation representative of the bulk
567 < solvent structural properties. Belch, et al \cite{Belch1985} simulated
568 < clusters of TIPS2 water surrounded by a hydrophobic bounding
569 < potential. The spherical hydrophobic boundary induced dangling
570 < hydrogen bonds at the surface that propagated deep into the cluster,
571 < affecting 70\% of the 100 molecules in the simulation. This result
572 < 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.
581 <
582 < In contrast, the Langevin Hull does not require that the orientation
583 < of molecules be fixed, nor does it utilize an explicitly hydrophobic
584 < boundary, orientational constraint or radial constraint. The number
585 < 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.
562 > As described above, the Langevin Hull does not require that the
563 > orientation of molecules be fixed, nor does it utilize an explicitly
564 > hydrophobic boundary, orientational constraint or radial constraint.
565 > Therefore, the orientational correlations of the molecules in a water
566 > cluster are of particular interest in testing this method.  Ideally,
567 > the water molecules on the surface of the cluster will have enough
568 > mobility into and out of the center of the cluster to maintain a
569 > bulk-like orientational distribution in the absence of orientational
570 > and radial constraints.  However, since the number of hydrogen bonding
571 > partners available to molecules on the exterior are limited, it is
572 > likely that there will be some effective hydrophobicity of the hull.
573  
574 < The orientation of a water molecule is described by
575 <
574 > To determine the extent of these effects demonstrated by the Langevin
575 > Hull, we examined the orientationations exhibited by SPC/E water in a
576 > cluster of 1372 molecules at 300 K and at pressures ranging from 1 -
577 > 1000 atm.  The orientational angle of a water molecule is described
578   \begin{equation}
579   \cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
580   \end{equation}
581 <
582 < 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}.
583 <
581 > where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of
582 > mass and the cluster center of mass and $\vec{\mu}_{i}$ is the vector
583 > bisecting the H-O-H angle of molecule {\it i} (See
584 > Fig. \ref{fig:coords}).
585   \begin{figure}
586   \includegraphics[width=\linewidth]{g_r_theta}
587 < \caption{Definition of coordinates}
588 < \label{coords}
587 > \caption{Orientation angle of the water molecules relative to the
588 >  center of the cluster.  Bulk-like distributions will result in
589 >  $\langle \cos \theta \rangle$ values close to zero.  If the hull
590 >  exhibits an overabundance of externally-oriented oxygen sites the
591 >  average orientation will be negative, while dangling hydrogen sites
592 >  will result in positive average orientations.}
593 > \label{fig:coords}
594   \end{figure}
595  
596 < 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).
597 <
596 > Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
597 > for molecules in the interior of the cluster (squares) and for
598 > molecules included in the convex hull (circles).
599   \begin{figure}
600   \includegraphics[width=\linewidth]{pAngle}
601 < \caption{SPC/E water clusters: only minor dewetting at the boundary}
602 < \label{pAngle}
601 > \caption{Distribution of $\cos{\theta}$ values for molecules on the
602 >  interior of the cluster (squares) and for those participating in the
603 >  convex hull (circles) at a variety of pressures.  The Langevin hull
604 >  exhibits minor dewetting behavior with exposed oxygen sites on the
605 >  hull water molecules.  The orientational preference for exposed
606 >  oxygen appears to be independent of applied pressure. }
607 > \label{fig:pAngle}
608   \end{figure}
609  
610 < 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.
610 > As expected, interior molecules (those not included in the convex
611 > hull) maintain a bulk-like structure with a uniform distribution of
612 > orientations. Molecules included in the convex hull show a slight
613 > preference for values of $\cos{\theta} < 0.$ These values correspond
614 > to molecules with oxygen directed toward the exterior of the cluster,
615 > forming a dangling hydrogen bond acceptor site.
616  
617 < 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.
617 > In the absence of an electrostatic contribution from the exterior
618 > bath, the orientational distribution of water molecules included in
619 > the Langevin Hull will slightly resemble the distribution at a neat
620 > water liquid/vapor interface.  Previous molecular dynamics simulations
621 > of SPC/E water \cite{Taylor1996} have shown that molecules at the
622 > liquid/vapor interface favor an orientation where one hydrogen
623 > protrudes from the liquid phase. This behavior is demonstrated by
624 > experiments \cite{Du1994} \cite{Scatena2001} showing that
625 > approximately one-quarter of water molecules at the liquid/vapor
626 > interface form dangling hydrogen bonds. The negligible preference
627 > shown in these cluster simulations could be removed through the
628 > introduction of an implicit solvent model, which would provide the
629 > missing electrostatic interactions between the cluster molecules and
630 > the surrounding temperature/pressure bath.
631  
632 < 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.
632 > The orientational preference exhibited by hull molecules in the
633 > Langevin hull is significantly weaker than the preference caused by an
634 > explicit hydrophobic bounding potential.  Additionally, the Langevin
635 > Hull does not require that the orientation of any molecules be fixed
636 > in order to maintain bulk-like structure, even at the cluster surface.
637  
638   \subsection{Heterogeneous nanoparticle / water mixtures}
639  
640   \section{Discussion}
641   \label{sec:discussion}
642  
643 + The Langevin Hull samples the isobaric-isothermal ensemble for
644 + non-periodic systems by coupling the system to an bath characterized
645 + by pressure, temperature, and solvent viscosity.  This enables the
646 + study of heterogeneous systems composed of materials of significantly
647 + different compressibilities.  Because the boundary is dynamically
648 + determined during the simulation and the molecules interacting with
649 + the boundary can change, the method and has minimal perturbations on
650 + the behavior of molecules at the edges of the simulation.  Further
651 + work on this method will involve implicit electrostatics at the
652 + boundary (which is missing in the current implementation) as well as
653 + more sophisticated treatments of the surface geometry (alpha
654 + shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight
655 + Cocone\cite{Dey:2003ts}). The non-convex hull geometries are
656 + significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull
657 + ($\mathcal{O}(N \log N)$), but would enable the use of hull volumes
658 + directly in computing the compressibility of the sample.
659 +
660   \section*{Appendix A: Computing Convex Hulls on Parallel Computers}
661  
662 + In order to use the Langevin Hull for simulations on parallel
663 + computers, one of the more difficult tasks is to compute the bounding
664 + surface, facets, and resistance tensors when the processors have
665 + incomplete information about the entire system's topology.  Most
666 + parallel decomposition methods assign primary responsibility for the
667 + motion of an atomic site to a single processor, and we can exploit
668 + this to efficiently compute the convex hull for the entire system.
669 +
670 + The basic idea involves splitting the point cloud into
671 + spatially-overlapping subsets and computing the convex hulls for each
672 + of the subsets.  The points on the convex hull of the entire system
673 + are all present on at least one of the subset hulls. The algorithm
674 + works as follows:
675 + \begin{enumerate}
676 + \item Each processor computes the convex hull for its own atomic sites
677 +  (left panel in Fig. \ref{fig:parallel}).
678 + \item The Hull vertices from each processor are passed out to all of
679 +  the processors, and each processor assembles a complete list of hull
680 +  sites (this is much smaller than the original number of points in
681 +  the point cloud).
682 + \item Each processor computes the global convex hull (right panel in
683 +  Fig. \ref{fig:parallel}) using only those points that are the union
684 +  of sites gathered from all of the subset hulls.  Delaunay
685 +  triangulation is then done to obtain the facets of the global hull.
686 + \end{enumerate}
687 +
688 + \begin{figure}
689 + \begin{centering}
690 + \includegraphics[width=\linewidth]{parallel}
691 + \caption{When the sites are distributed among many nodes for parallel
692 +  computation, the processors first compute the convex hulls for their
693 +  own sites (dashed lines in left panel). The positions of the sites
694 +  that make up the convex hulls are then communicated to all
695 +  processors (middle panel).  The convex hull of the system (solid line in right panel) is the convex hull of the points on the hulls for all
696 +  processors.}
697 + \label{fig:parallel}
698 + \end{centering}
699 + \label{fig:parallel}
700 + \end{figure}
701 +
702 + The individual hull operations scale with
703 + $\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total
704 + number of sites, and $p$ is the number of processors.  These local
705 + hull operations create a set of $p$ hulls each with approximately
706 + $\frac{n}{3pr}$ sites (for a cluster of radius $r$). The worst-case
707 + communication cost for using a ``gather'' operation to distribute this
708 + information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n
709 +  \beta (p-1)}{3 r p^2})$, while the final computation of the system
710 + hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$.
711 +
712 + For a large number of atoms on a moderately parallel machine, the
713 + total costs are dominated by the computations of the individual hulls,
714 + and communication of these hulls to so the Langevin hull sees roughly
715 + linear speed-up with increasing processor counts.
716 +
717   \section*{Acknowledgments}
718   Support for this project was provided by the
719   National Science Foundation under grant CHE-0848243. Computational

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines