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 3688 by gezelter, Thu Nov 11 17:13:49 2010 UTC vs.
Revision 3689 by gezelter, Thu Nov 11 22:17:37 2010 UTC

# Line 117 | Line 117 | expensive.  A 62 $\AA^3$ box of water solvating a mode
117   higher than desirable, or unreasonable system sizes to avoid this
118   effect.  For example, calculations using typical hydration boxes
119   solvating a protein under periodic boundary conditions are quite
120 < expensive.  A 62 $\AA^3$ box of water solvating a moderately small
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 protein concentrations in the cell are on the order of
125 < 160-310 mg/ml,\cite{Brown1991195} and the factor of 20 difference
126 < between simulations and the cellular environment may have significant
127 < effects on the structure and dynamics of simulated protein structures.
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  
129
131   \subsection*{Boundary Methods}
132   There have been a number of approaches to handle simulations of
133   explicitly non-periodic systems that focus on constant or
# Line 138 | Line 139 | used widely for protein simulations. [CITATIONS NEEDED
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. [CITATIONS NEEDED]
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 149 | 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 180 | 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 196 | 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 221 | Line 222 | The Langevin Hull uses an external bath at a fixed con
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 300 | 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 330 | 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 454 | 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{Compressibility of gold nanoparticles}
459 > \subsection{Bulk Modulus of gold nanoparticles}
460  
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
464 < and related multi-body force fields.  In particular, the quantum
465 < Sutton-Chen potential gets nearly quantitative agreement with the
466 < experimental bulk modulus values, and makes a good first test of how
467 < the Langevin Hull will perform at large applied pressures.
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 that closely
473 < resembles the Lennard Jones potential,
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   \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] ,
# Line 490 | Line 492 | metals.\cite{PhysRevB.59.3527}
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}
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 have yielded values of
500 < 175.53 GPa.\cite{XXX} Using the same force field, we have performed a
501 < series of relatively short (200 ps) simulations on 40 \r{A} radius
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 GPa to XXX.  We obtain a value of 177.547 GPa
504 < for the bulk modulus for gold using this echnique.
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]{stacked}
# Line 508 | Line 512 | for the bulk modulus for gold using this echnique.
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{pressureResponse}
515 > \label{fig:pressureResponse}
516   \end{figure}
517  
518 < \begin{equation}
519 < \kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial
520 <    P}\right)
521 < \end{equation}
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 526 | 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.
529 Compressibility values from all references are for applied pressures
530 within the range 1 - 1000 atm.
539  
540   \begin{figure}
541   \includegraphics[width=\linewidth]{new_isothermalN}
# Line 537 | 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 553 | Line 561 | geometries which include large volumes of empty space.
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 574 | 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. However, 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.  Early non-periodic
607 < simulation methods (e.g. hydrophobic boundary potentials) induced
608 < spurious orientational correlations deep within the simulated
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
# Line 602 | Line 613 | hydrophobic boundary, orientational constraint or radi
613  
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, orientational constraint or radial constraint.
617 < Therefore, the orientational correlations of the molecules in a water
618 < cluster are of particular interest in testing this method.  Ideally,
619 < the water molecules on the surface of the cluster will have enough
620 < mobility into and out of the center of the cluster to maintain a
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 some effective hydrophobicity of the hull.
624 > likely that there will be an effective hydrophobicity of the hull.
625  
626 < To determine the extent of these effects demonstrated by the Langevin
627 < Hull, we examined the orientationations exhibited by SPC/E water in a
628 < cluster of 1372 molecules at 300 K and at pressures ranging from 1 -
629 < 1000 atm.  The orientational angle of a water molecule is described
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 distributions
636 < will result in $\langle \cos \theta \rangle$ values close to zero.  If
637 < the hull exhibits an overabundance of externally-oriented oxygen sites
638 < the average orientation will be negative, while dangling hydrogen
639 < sites will result in positive average orientations.
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   Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values
643   for molecules in the interior of the cluster (squares) and for
# Line 671 | Line 683 | in order to maintain bulk-like structure, even at the
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 an bath characterized
717 < by pressure, temperature, and solvent viscosity.  This enables the
718 < study of heterogeneous systems composed of materials of significantly
719 < different compressibilities.  Because the boundary is dynamically
720 < determined during the simulation and the molecules interacting with
721 < the boundary can change, the method and has minimal perturbations on
722 < the behavior of molecules at the edges of the simulation.  Further
723 < work on this method will involve implicit electrostatics at the
724 < boundary (which is missing in the current implementation) as well as
725 < more sophisticated treatments of the surface geometry (alpha
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
# Line 695 | Line 734 | surface, facets, and resistance tensors when the proce
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 processors have
738 < incomplete information about the entire system's topology.  Most
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.

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines