17 |
|
\setlength{\abovecaptionskip}{20 pt} |
18 |
|
\setlength{\belowcaptionskip}{30 pt} |
19 |
|
|
20 |
< |
\bibpunct{[}{]}{,}{s}{}{;} |
20 |
> |
\bibpunct{}{}{,}{s}{}{;} |
21 |
|
\bibliographystyle{achemso} |
22 |
|
|
23 |
|
\begin{document} |
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 |
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 |
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. |
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 |
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 |
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 |
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 |
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). |
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 |
|
|
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 |
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}, |
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 |
|
|
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 |
|
|
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 |
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 |
< |
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. |
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 |
< |
Previous molecular dynamics simulations |
704 |
< |
of SPC/E water using 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 while we are unable to comment on the orientational preference of vapor phase molecules in a Langevin Hull simulation, we achieve good agreement for the orientation of liquid phase molecules at the interface. |
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 and |
726 |
< |
the density of gold at the surface as a function of pressure. |
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 |
< |
Higher applied pressures de-structure the outermost layer of the gold nanoparticle and the water at the metal/water interface. Increased pressure shows more overlap of the gold and water densities, indicating a less well-defined interfacial surface. |
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 |
|
\begin{figure} |
739 |
|
\includegraphics[width=\linewidth]{RhoR} |
740 |
< |
\caption{Densities of gold and water at the nanoparticle surface. Higher applied pressures de-structure both the gold nanoparticle surface and water at the metal/water interface.} |
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} |