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 3693 by kstocke1, Mon Nov 15 16:46:53 2010 UTC vs.
Revision 3723 by gezelter, Fri Feb 18 14:08:15 2011 UTC

# Line 17 | Line 17
17   \setlength{\abovecaptionskip}{20 pt}
18   \setlength{\belowcaptionskip}{30 pt}
19  
20 < \bibpunct{[}{]}{,}{s}{}{;}
20 > \bibpunct{}{}{,}{s}{}{;}
21   \bibliographystyle{achemso}
22  
23   \begin{document}
# Line 42 | Line 42 | Notre Dame, Indiana 46556}
42    hull surrounding the system.  A Langevin thermostat is also applied
43    to the facets to mimic contact with an external heat bath. This new
44    method, the ``Langevin Hull'', can handle heterogeneous mixtures of
45 <  materials with different compressibilities.  These are systems that
46 <  are problematic for traditional affine transform methods.  The
47 <  Langevin Hull does not suffer from the edge effects of boundary
48 <  potential methods, and allows realistic treatment of both external
49 <  pressure and thermal conductivity due to the presence of an implicit
50 <  solvent.  We apply this method to several different systems
51 <  including bare metal nanoparticles, nanoparticles in an explicit
52 <  solvent, as well as clusters of liquid water. The predicted
53 <  mechanical properties of these systems are in good agreement with
54 <  experimental data and previous simulation work.
45 >  materials with different compressibilities.  These systems are
46 >  problematic for traditional affine transform methods.  The Langevin
47 >  Hull does not suffer from the edge effects of boundary potential
48 >  methods, and allows realistic treatment of both external pressure
49 >  and thermal conductivity due to the presence of an implicit solvent.
50 >  We apply this method to several different systems including bare
51 >  metal nanoparticles, nanoparticles in an explicit solvent, as well
52 >  as clusters of liquid water. The predicted mechanical properties of
53 >  these systems are in good agreement with experimental data and
54 >  previous simulation work.
55   \end{abstract}
56  
57   \newpage
# Line 121 | Line 121 | effective protein concentration of 100 mg/mL.\cite{Ast
121   protein like hen egg white lysozyme (PDB code: 1LYZ) yields an
122   effective protein concentration of 100 mg/mL.\cite{Asthagiri20053300}
123  
124 < {\it Yotal} protein concentrations in the cell are typically on the
124 > {\it Total} protein concentrations in the cell are typically on the
125   order of 160-310 mg/ml,\cite{Brown1991195} and individual proteins
126   have concentrations orders of magnitude lower than this in the
127   cellular environment. The effective concentrations of single proteins
128   in simulations may have significant effects on the structure and
129 < dynamics of simulated structures.
129 > dynamics of simulated systems.
130  
131   \subsection*{Boundary Methods}
132   There have been a number of approaches to handle simulations of
# Line 251 | Line 251 | simulation.
251   \caption{The external temperature and pressure bath interacts only
252    with those atoms on the convex hull (grey surface).  The hull is
253    computed dynamically at each time step, and molecules can move
254 <  between the interior (Newtonian) region and the Langevin hull.}
254 >  between the interior (Newtonian) region and the Langevin Hull.}
255   \label{fig:hullSample}
256   \end{figure}
257  
# Line 266 | Line 266 | equation of motion is modified with an external force,
266   potential energy.  For atoms on the exterior of the cluster
267   (i.e. those that occupy one of the vertices of the convex hull), the
268   equation of motion is modified with an external force, ${\mathbf
269 <  F}_i^{\mathrm ext}$,
269 >  F}_i^{\mathrm ext}$:
270   \begin{equation}
271   m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}.
272   \end{equation}
# Line 378 | Line 378 | integrator in our code, OpenMD.\cite{Meineke2005,openm
378   configurations, so this appears to be a reasonably good approximation.
379  
380   We have implemented this method by extending the Langevin dynamics
381 < integrator in our code, OpenMD.\cite{Meineke2005,openmd}  At each
381 > integrator in our code, OpenMD.\cite{Meineke2005,open_md}  At each
382   molecular dynamics time step, the following process is carried out:
383   \begin{enumerate}
384   \item The standard inter-atomic forces ($\nabla_iU$) are computed.
# Line 400 | Line 400 | using calls to the qhull library.\cite{Qhull} There is
400   \item Atomic positions and velocities are propagated.
401   \end{enumerate}
402   The Delaunay triangulation and computation of the convex hull are done
403 < using calls to the qhull library.\cite{Qhull} There is a minimal
403 > using calls to the qhull library.\cite{Q_hull} There is a minimal
404   penalty for computing the convex hull and resistance tensors at each
405   step in the molecular dynamics simulation (roughly 0.02 $\times$ cost
406   of a single force evaluation), and the convex hull is remarkably easy
# Line 412 | Line 412 | heterogeneous mixture (gold nanoparticles in a water d
412   To test the new method, we have carried out simulations using the
413   Langevin Hull on: 1) a crystalline system (gold nanoparticles), 2) a
414   liquid droplet (SPC/E water),\cite{Berendsen1987} and 3) a
415 < heterogeneous mixture (gold nanoparticles in a water droplet). In each
416 < case, we have computed properties that depend on the external applied
417 < pressure.  Of particular interest for the single-phase systems is the
418 < isothermal compressibility,
415 > heterogeneous mixture (gold nanoparticles in an SPC/E water droplet). In each case, we have computed properties that depend on the external applied pressure. Of particular interest for the single-phase systems is the isothermal compressibility,
416   \begin{equation}
417   \kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right
418   )_{T}.
# Line 424 | Line 421 | is not well-defined.  In order to compute the compress
421  
422   One problem with eliminating periodic boundary conditions and
423   simulation boxes is that the volume of a three-dimensional point cloud
424 < is not well-defined.  In order to compute the compressibility of a
424 > is not well-defined. In order to compute the compressibility of a
425   bulk material, we make an assumption that the number density, $\rho =
426 < \frac{N}{V}$, is uniform within some region of the point cloud.  The
426 > \frac{N}{V}$, is uniform within some region of the point cloud. The
427   compressibility can then be expressed in terms of the average number
428   of particles in that region,
429   \begin{equation}
430   \kappa_{T} = -\frac{1}{N} \left ( \frac{\partial N}{\partial P} \right
431 < )_{T}
431 > )_{T}.
432   \label{eq:BMN}
433   \end{equation}
434 < The region we used is a spherical volume of 10 \AA\ radius centered in
435 < the middle of the cluster. $N$ is the average number of molecules
434 > The region we used is a spherical volume of 20 \AA\ radius centered in
435 > the middle of the cluster with a roughly 25 \AA\ radius. $N$ is the average number of molecules
436   found within this region throughout a given simulation. The geometry
437 < and size of the region is arbitrary, and any bulk-like portion of the
438 < cluster can be used to compute the compressibility.
437 > of the region is arbitrary, and any bulk-like portion of the
438 > cluster can be used to compute the compressibility.
439  
440   One might assume that the volume of the convex hull could simply be
441   taken as the system volume $V$ in the compressibility expression
# Line 484 | Line 481 | pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate o
481   interactions between the pseudoatom cores. The $\sqrt{\rho_i}$ term in
482   Eq. (\ref{eq:SCP1}) is an attractive many-body potential that models
483   the interactions between the valence electrons and the cores of the
484 < pseudo-atoms. $D_{ij}$, $D_{ii}$ set the appropriate overall energy
484 > pseudo-atoms. $D_{ij}$ and $D_{ii}$ set the appropriate overall energy
485   scale, $c_i$ scales the attractive portion of the potential relative
486   to the repulsive interaction and $\alpha_{ij}$ is a length parameter
487   that assures a dimensionless form for $\rho$. These parameters are
# Line 492 | Line 489 | metals.\cite{PhysRevB.59.3527,QSC}
489   energy, and elastic moduli for FCC transition metals. The quantum
490   Sutton-Chen (QSC) formulation matches these properties while including
491   zero-point quantum corrections for different transition
492 < metals.\cite{PhysRevB.59.3527,QSC}
492 > metals.\cite{PhysRevB.59.3527,QSC2}
493  
494   In bulk gold, the experimentally-measured value for the bulk modulus
495   is 180.32 GPa, while previous calculations on the QSC potential in
496   periodic-boundary simulations of the bulk crystal have yielded values
497 < of 175.53 GPa.\cite{QSC} Using the same force field, we have performed
498 < a series of 1 ns simulations on 40 \AA~ radius
499 < nanoparticles under the Langevin Hull at a variety of applied
500 < pressures ranging from 0 -- 10 GPa.  We obtain a value of 177.55 GPa
501 < for the bulk modulus of gold using this technique, in close agreement
502 < with both previous simulations and the experimental bulk modulus of
503 < gold.
497 > of 175.53 GPa.\cite{QSC2} Using the same force field, we have
498 > performed a series of 1 ns simulations on gold nanoparticles of three
499 > different radii: 20 \AA~ (1985 atoms), 30 \AA~ (6699 atoms), and 40
500 > \AA~ (15707 atoms) utilizing the Langevin Hull at a variety of applied
501 > pressures ranging from 0 -- 10 GPa.  For the 40 \AA~ radius
502 > nanoparticle we obtain a value of 177.55 GPa for the bulk modulus of
503 > gold, in close agreement with both previous simulations and the
504 > experimental bulk modulus reported for gold single
505 > crystals.\cite{Collard1991} The smaller gold nanoparticles (30 and 20
506 > \AA~ radii) have calculated bulk moduli of 215.58 and 208.86 GPa,
507 > respectively, indicating that smaller nanoparticles are somewhat
508 > stiffer (less compressible) than the larger nanoparticles.  This
509 > stiffening of the small nanoparticles may be related to their high
510 > degree of surface curvature, resulting in a lower coordination number
511 > of surface atoms relative to the the surface atoms in the 40 \AA~
512 > radius particle.
513  
514 + We obtain a gold lattice constant of 4.051 \AA~ using the Langevin
515 + Hull at 1 atm, close to the experimentally-determined value for bulk
516 + gold and the value for gold simulated using the QSC potential and
517 + periodic boundary conditions (4.079 \AA~ and 4.088\AA~,
518 + respectively).\cite{QSC2} The slightly smaller calculated lattice
519 + constant is most likely due to the presence of surface tension in the
520 + non-periodic Langevin Hull cluster, an effect absent from a bulk
521 + simulation. The specific heat of a 40 \AA~ gold nanoparticle under the
522 + Langevin Hull at 1 atm is 24.914 $\mathrm {\frac{J}{mol \, K}}$, which
523 + compares very well with the experimental value of 25.42 $\mathrm
524 + {\frac{J}{mol \, K}}$.
525 +
526   \begin{figure}
527   \includegraphics[width=\linewidth]{stacked}
528   \caption{The response of the internal pressure and temperature of gold
# Line 521 | Line 539 | the total surface area of the cluter exposed to the ba
539   temperature respond to the Langevin Hull for nanoparticles that were
540   initialized far from the target pressure and temperature.  As
541   expected, the rate at which thermal equilibrium is achieved depends on
542 < the total surface area of the cluter exposed to the bath as well as
542 > the total surface area of the cluster exposed to the bath as well as
543   the bath viscosity.  Pressure that is applied suddenly to a cluster
544   can excite breathing vibrations, but these rapidly damp out (on time
545 < scales of 30-50 ps).
545 > scales of 30 -- 50 ps).
546  
547   \subsection{Compressibility of SPC/E water clusters}
548  
# Line 533 | Line 551 | Langevin Hull simulations for pressures between 1 and
551   ensembles) have yielded values for the isothermal compressibility that
552   agree well with experiment.\cite{Fine1973} The results of two
553   different approaches for computing the isothermal compressibility from
554 < Langevin Hull simulations for pressures between 1 and 6500 atm are
554 > Langevin Hull simulations for pressures between 1 and 3000 atm are
555   shown in Fig. \ref{fig:compWater} along with compressibility values
556   obtained from both other SPC/E simulations and experiment.
557  
# Line 548 | Line 566 | pressures.  The reason for this deviation is quite sim
566   and previous simulation work throughout the 1 -- 1000 atm pressure
567   regime.  Compressibilities computed using the Hull volume, however,
568   deviate dramatically from the experimental values at low applied
569 < pressures.  The reason for this deviation is quite simple; at low
569 > pressures.  The reason for this deviation is quite simple: at low
570   applied pressures, the liquid is in equilibrium with a vapor phase,
571   and it is entirely possible for one (or a few) molecules to drift away
572   from the liquid cluster (see Fig. \ref{fig:coneOfShame}).  At low
# Line 578 | Line 596 | volume,\cite{Debenedetti1986},
596   different pressures must be done to compute the first derivatives.  It
597   is also possible to compute the compressibility using the fluctuation
598   dissipation theorem using either fluctuations in the
599 < volume,\cite{Debenedetti1986},
599 > volume,\cite{Debenedetti1986}
600   \begin{equation}
601   \kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle
602      V \right \rangle ^{2}}{V \, k_{B} \, T},
# Line 588 | Line 606 | fixed region,
606   fixed region,
607   \begin{equation}
608   \kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle
609 <    N \right \rangle ^{2}}{N \, k_{B} \, T},
609 >    N \right \rangle ^{2}}{N \, k_{B} \, T}.
610   \label{eq:BMNfluct}
611   \end{equation}
612   Thus, the compressibility of each simulation can be calculated
# Line 597 | Line 615 | compressibilities.
615   effects of the empty space due to the vapor phase; for this reason, we
616   recommend using the number density (Eq. \ref{eq:BMN}) or number
617   density fluctuations (Eq. \ref{eq:BMNfluct}) for computing
618 < compressibilities.
618 > compressibilities.  We obtained the results in
619 > Fig. \ref{fig:compWater} using a sampling radius that was
620 > approximately 80\% of the mean distance between the center of mass of
621 > the cluster and the hull atoms.  This ratio of sampling radius to
622 > average hull radius excludes the problematic vapor phase on the
623 > outside of the cluster while including enough of the liquid phase to
624 > avoid poor statistics due to fluctuating local densities.
625  
626 + A comparison of the oxygen-oxygen radial distribution functions for
627 + SPC/E water simulated using both the Langevin Hull and more
628 + traditional periodic boundary methods -- both at 1 atm and 300K --
629 + reveals an understructuring of water in the Langevin Hull that
630 + manifests as a slight broadening of the solvation shells.  This effect
631 + may be due to the introduction of a liquid-vapor interface in the
632 + Langevin Hull simulations (an interface which is missing in most
633 + periodic simulations of bulk water).  Vapor-phase molecules contribute
634 + a small but nearly flat portion of the radial distribution function.
635 +
636   \subsection{Molecular orientation distribution at cluster boundary}
637  
638   In order for a non-periodic boundary method to be widely applicable,
# Line 607 | Line 641 | fixing and characterizing the effects of artifical bou
641   methods (e.g. hydrophobic boundary potentials) induced spurious
642   orientational correlations deep within the simulated
643   system.\cite{Lee1984,Belch1985} This behavior spawned many methods for
644 < fixing and characterizing the effects of artifical boundaries
644 > fixing and characterizing the effects of artificial boundaries
645   including methods which fix the orientations of a set of edge
646   molecules.\cite{Warshel1978,King1989}
647  
# Line 616 | Line 650 | the water molecules on the surfaces of the clusterss w
650   hydrophobic boundary, or orientational or radial constraints.
651   Therefore, the orientational correlations of the molecules in water
652   clusters are of particular interest in testing this method.  Ideally,
653 < the water molecules on the surfaces of the clusterss will have enough
653 > the water molecules on the surfaces of the clusters will have enough
654   mobility into and out of the center of the cluster to maintain
655   bulk-like orientational distribution in the absence of orientational
656   and radial constraints.  However, since the number of hydrogen bonding
# Line 646 | Line 680 | molecules included in the convex hull (circles).
680   \includegraphics[width=\linewidth]{pAngle}
681   \caption{Distribution of $\cos{\theta}$ values for molecules on the
682    interior of the cluster (squares) and for those participating in the
683 <  convex hull (circles) at a variety of pressures.  The Langevin hull
683 >  convex hull (circles) at a variety of pressures.  The Langevin Hull
684    exhibits minor dewetting behavior with exposed oxygen sites on the
685    hull water molecules.  The orientational preference for exposed
686    oxygen appears to be independent of applied pressure. }
# Line 658 | Line 692 | forming a dangling hydrogen bond acceptor site.
692   orientations. Molecules included in the convex hull show a slight
693   preference for values of $\cos{\theta} < 0.$ These values correspond
694   to molecules with oxygen directed toward the exterior of the cluster,
695 < forming a dangling hydrogen bond acceptor site.
695 > forming dangling hydrogen bond acceptor sites.
696  
697 < Previous molecular dynamics simulations
698 < of SPC/E water that use periodic boundary conditions have shown that molecules on the liquid side of the liquid/vapor interface favor a similar orientation where oxygen is directed away from the bulk.\cite{Taylor1996} These simulations had both a liquid phase and a well-defined vapor phase in equilibrium and showed that vapor molecules generally had one hydrogen protruding from the surface, forming a dangling hydrogen bond donor. Our water cluster simulations do not have a true lasting vapor phase, but rather a few transient molecules that leave the liquid droplet. Thus we are unable to comment on the orientational preference of vapor phase molecules in a Langevin Hull simulation.
697 > The orientational preference exhibited by water molecules on the hull
698 > is significantly weaker than the preference caused by an explicit
699 > hydrophobic bounding potential.  Additionally, the Langevin Hull does
700 > not require that the orientation of any molecules be fixed in order to
701 > maintain bulk-like structure, even near the cluster surface.
702  
703 < However, the orientational preference exhibited by liquid phase hull molecules in the Langevin hull is significantly weaker than the preference caused by an explicit hydrophobic bounding potential.  Additionally, the Langevin Hull does not require that the orientation of any molecules be fixed in order to maintain bulk-like structure, even at the cluster surface.
703 > Previous molecular dynamics simulations of SPC/E liquid / vapor
704 > interfaces using periodic boundary conditions have shown that
705 > molecules on the liquid side of interface favor a similar orientation
706 > where oxygen is directed away from the bulk.\cite{Taylor1996} These
707 > simulations had well-defined liquid and vapor phase regions
708 > equilibrium and it was observed that {\it vapor} molecules generally
709 > had one hydrogen protruding from the surface, forming a dangling
710 > hydrogen bond donor. Our water clusters do not have a true vapor
711 > region, but rather a few transient molecules that leave the liquid
712 > droplet (and which return to the droplet relatively quickly).
713 > Although we cannot obtain an orientational preference of vapor phase
714 > molecules in a Langevin Hull simulation, but we do agree with previous
715 > estimates of the orientation of {\it liquid phase} molecules at the
716 > interface.
717  
718   \subsection{Heterogeneous nanoparticle / water mixtures}
719  
720 < To further test the method, we simulated gold nanopartices ($r = 18$
721 < \AA) solvated by explicit SPC/E water clusters using the Langevin
722 < hull.  This was done at pressures of 1, 2, 5, 10, 20, 50, 100 and 200 atm
723 < in order to observe the effects of pressure on the ordering of water
724 < ordering at the surface.  In Fig. \ref{fig:RhoR} we show the density
725 < of water adjacent to the surface as a function of pressure, as well as
726 < the orientational ordering of water at the surface of the
727 < nanoparticle.
720 > To further test the method, we simulated gold nanoparticles ($r = 18$
721 > \AA~, 1433 atoms) solvated by explicit SPC/E water clusters (5000
722 > molecules) using a model for the gold / water interactions that has
723 > been used by Dou {\it et. al.} for investigating the separation of
724 > water films near hot metal surfaces.\cite{ISI:000167766600035} The
725 > Langevin Hull was used to sample pressures of 1, 2, 5, 10, 20, 50, 100
726 > and 200 atm, while all simulations were done at a temperature of 300
727 > K.  At these temperatures and pressures, there is no observed
728 > separation of the water film from the surface.
729  
730 < \begin{figure}
730 > In Fig. \ref{fig:RhoR} we show the density of water and gold as a
731 > function of the distance from the center of the nanoparticle.  Higher
732 > applied pressures appear to destroy structural correlations in the
733 > outermost monolayer of the gold nanoparticle as well as in the water
734 > at the near the metal / water interface.  Simulations at increased
735 > pressures exhibit significant overlap of the gold and water densities,
736 > indicating a less well-defined interfacial surface.
737  
738 < \caption{interesting plot showing cluster behavior}
738 > \begin{figure}
739 > \includegraphics[width=\linewidth]{RhoR}
740 > \caption{Density profiles of gold and water at the nanoparticle
741 >  surface. Each curve has been normalized by the average density in
742 >  the bulk-like region available to the corresponding material.
743 >  Higher applied pressures de-structure both the gold nanoparticle
744 >  surface and water at the metal/water interface.}
745   \label{fig:RhoR}
746   \end{figure}
747  
748 < At higher pressures, problems with the gold - water interaction
749 < potential became apparent.  The model we are using (due to Spohr) was
750 < intended for relatively low pressures; it utilizes both shifted Morse
751 < and repulsive Morse potentials to model the Au/O and Au/H
752 < interactions, respectively.  The repulsive wall of the Morse potential
753 < does not diverge quickly enough at short distances to prevent water
754 < from diffusing into the center of the gold nanoparticles.  This
755 < behavior is likely not a realistic description of the real physics of
756 < the situation.  A better model of the gold-water adsorption behavior
757 < appears to require harder repulsive walls to prevent this behavior.
748 > At even higher pressures (500 atm and above), problems with the metal
749 > - water interaction potential became quite clear.  The model we are
750 > using appears to have been parameterized for relatively low pressures;
751 > it utilizes both shifted Morse and repulsive Morse potentials to model
752 > the Au/O and Au/H interactions, respectively.  The repulsive wall of
753 > the Morse potential does not diverge quickly enough at short distances
754 > to prevent water from diffusing into the center of the gold
755 > nanoparticles.  This behavior is likely not a realistic description of
756 > the real physics of the situation.  A better model of the gold-water
757 > adsorption behavior would require harder repulsive walls to prevent
758 > this behavior.
759  
760   \section{Discussion}
761   \label{sec:discussion}
# Line 766 | Line 830 | and communication of these hulls to create the Langevi
830  
831   For a large number of atoms on a moderately parallel machine, the
832   total costs are dominated by the computations of the individual hulls,
833 < and communication of these hulls to create the Langevin hull sees roughly
833 > and communication of these hulls to create the Langevin Hull sees roughly
834   linear speed-up with increasing processor counts.
835  
836   \section*{Acknowledgments}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines