| 39 |
|
\begin{abstract} |
| 40 |
|
We have developed a new isobaric-isothermal (NPT) algorithm which |
| 41 |
|
applies an external pressure to the facets comprising the convex |
| 42 |
< |
hull surrounding the objects in the system. Additionally, a Langevin |
| 43 |
< |
thermostat is applied to facets of the hull to mimic contact with an |
| 44 |
< |
external heat bath. This new method, the ``Langevin Hull'', performs |
| 45 |
< |
better than traditional affine transform methods for systems |
| 46 |
< |
containing heterogeneous mixtures of materials with different |
| 42 |
> |
hull surrounding the system. A Langevin thermostat is also applied |
| 43 |
> |
to facets of the hull to mimic contact with an external heat |
| 44 |
> |
bath. This new method, the ``Langevin Hull'', performs better than |
| 45 |
> |
traditional affine transform methods for systems containing |
| 46 |
> |
heterogeneous mixtures of materials with different |
| 47 |
|
compressibilities. It does not suffer from the edge effects of |
| 48 |
|
boundary potential methods, and allows realistic treatment of both |
| 49 |
|
external pressure and thermal conductivity to an implicit solvent. |
| 50 |
|
We apply this method to several different systems including bare |
| 51 |
< |
nanoparticles, nanoparticles in an explicit solvent, as well as |
| 52 |
< |
clusters of liquid water and ice. The predicted mechanical and |
| 53 |
< |
thermal properties of these systems are in good agreement with |
| 54 |
< |
experimental data. |
| 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 |
| 66 |
|
\section{Introduction} |
| 67 |
|
|
| 68 |
|
The most common molecular dynamics methods for sampling configurations |
| 69 |
< |
of an isobaric-isothermal (NPT) ensemble attempt to maintain a target |
| 70 |
< |
pressure in a simulation by coupling the volume of the system to an |
| 71 |
< |
extra degree of freedom, the {\it barostat}. These methods require |
| 72 |
< |
periodic boundary conditions, because when the instantaneous pressure |
| 73 |
< |
in the system differs from the target pressure, the volume is |
| 74 |
< |
typically reduced or expanded using {\it affine transforms} of the |
| 75 |
< |
system geometry. An affine transform scales both the box lengths as |
| 76 |
< |
well as the scaled particle positions (but not the sizes of the |
| 77 |
< |
particles). The most common constant pressure methods, including the |
| 78 |
< |
Melchionna modification\cite{Melchionna1993} to the |
| 79 |
< |
Nos\'e-Hoover-Andersen equations of |
| 80 |
< |
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen |
| 81 |
< |
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin |
| 82 |
< |
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize coordinate |
| 83 |
< |
transformation to adjust the box volume. As long as the material in |
| 84 |
< |
the simulation box is essentially a bulk-like liquid which has a |
| 85 |
< |
relatively uniform compressibility, the standard affine transform |
| 86 |
< |
approach provides an excellent way of adjusting the volume of the |
| 87 |
< |
system and applying pressure directly via the interactions between |
| 88 |
< |
atomic sites. |
| 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 |
| 83 |
> |
coordinate transformation to adjust the box volume. As long as the |
| 84 |
> |
material in the simulation box is essentially a bulk-like liquid which |
| 85 |
> |
has a relatively uniform compressibility, the standard affine |
| 86 |
> |
transform approach provides an excellent way of adjusting the volume |
| 87 |
> |
of the system and applying pressure directly via the interactions |
| 88 |
> |
between atomic sites. |
| 89 |
|
|
| 90 |
< |
The problem with this approach becomes apparent when the material |
| 91 |
< |
being simulated is an inhomogeneous mixture in which portions of the |
| 92 |
< |
simulation box are incompressible relative to other portions. |
| 93 |
< |
Examples include simulations of metallic nanoparticles in liquid |
| 94 |
< |
environments, proteins at interfaces, as well as other multi-phase or |
| 90 |
> |
One problem with this approach appears when the system being simulated |
| 91 |
> |
is an inhomogeneous mixture in which portions of the simulation box |
| 92 |
> |
are incompressible relative to other portions. Examples include |
| 93 |
> |
simulations of metallic nanoparticles in liquid environments, proteins |
| 94 |
> |
at ice / water interfaces, as well as other heterogeneous or |
| 95 |
|
interfacial environments. In these cases, the affine transform of |
| 96 |
|
atomic coordinates will either cause numerical instability when the |
| 97 |
< |
sites in the incompressible medium collide with each other, or lead to |
| 98 |
< |
inefficient sampling of system volumes if the barostat is set slow |
| 99 |
< |
enough to avoid the instabilities in the incompressible region. |
| 97 |
> |
sites in the incompressible medium collide with each other, or will |
| 98 |
> |
lead to inefficient sampling of system volumes if the barostat is set |
| 99 |
> |
slow enough to avoid the instabilities in the incompressible region. |
| 100 |
|
|
| 101 |
|
\begin{figure} |
| 102 |
|
\includegraphics[width=\linewidth]{AffineScale2} |
| 113 |
|
One may also wish to avoid affine transform periodic boundary methods |
| 114 |
|
to simulate {\it explicitly non-periodic systems} under constant |
| 115 |
|
pressure conditions. The use of periodic boxes to enforce a system |
| 116 |
< |
volume either requires effective solute concentrations that are much |
| 116 |
> |
volume requires either effective solute concentrations that are much |
| 117 |
|
higher than desirable, or unreasonable system sizes to avoid this |
| 118 |
< |
effect. For example, calculations using typical hydration shells |
| 118 |
> |
effect. For example, calculations using typical hydration shells |
| 119 |
|
solvating a protein under periodic boundary conditions are quite |
| 120 |
|
expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL |
| 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 |
| 134 |
|
simulations. [CITATIONS NEEDED] |
| 135 |
|
|
| 136 |
|
The electrostatic and dispersive behavior near the boundary has long |
| 137 |
< |
been a cause for concern. King and Warshel introduced a surface |
| 138 |
< |
constrained all-atom solvent (SCAAS) which included polarization |
| 139 |
< |
effects of a fixed spherical boundary to mimic bulk-like behavior |
| 140 |
< |
without periodic boundaries.\cite{king:3647} In the SCAAS model, a |
| 141 |
< |
layer of fixed solvent molecules surrounds the solute and any explicit |
| 142 |
< |
solvent, and this in turn is surrounded by a continuum dielectric. |
| 143 |
< |
MORE HERE. WHAT DID THEY FIND? |
| 144 |
< |
|
| 145 |
< |
Beglov and Roux developed a boundary model in which the hard sphere |
| 146 |
< |
boundary has a radius that varies with the instantaneous configuration |
| 147 |
< |
of the solute (and solvent) molecules.\cite{beglov:9050} This model |
| 148 |
< |
contains a clear pressure and surface tension contribution to the free |
| 149 |
< |
energy which XXX. |
| 137 |
> |
been a cause for concern when performing simulations of explicitly |
| 138 |
> |
non-periodic systems. Early work led to the surface constrained soft |
| 139 |
> |
sphere dipole model (SCSSD)\cite{Warshel1978} in which the surface |
| 140 |
> |
molecules are fixed in a random orientation representative of the bulk |
| 141 |
> |
solvent structural properties. Belch {\it et al.}\cite{Belch1985} |
| 142 |
> |
simulated clusters of TIPS2 water surrounded by a hydrophobic bounding |
| 143 |
> |
potential. The spherical hydrophobic boundary induced dangling |
| 144 |
> |
hydrogen bonds at the surface that propagated deep into the cluster, |
| 145 |
> |
affecting most of molecules in the simulation. This result echoes an |
| 146 |
> |
earlier study which showed that an extended planar hydrophobic surface |
| 147 |
> |
caused orientational preference at the surface which extended |
| 148 |
> |
relatively deep (7 \r{A}) into the liquid simulation |
| 149 |
> |
cell.\cite{Lee1984} The surface constrained all-atom solvent (SCAAS) |
| 150 |
> |
model \cite{King1989} improved upon its SCSSD predecessor. The SCAAS |
| 151 |
> |
model utilizes a polarization constraint which is applied to the |
| 152 |
> |
surface molecules to maintain bulk-like structure at the cluster |
| 153 |
> |
surface. A radial constraint is used to maintain the desired bulk |
| 154 |
> |
density of the liquid. Both constraint forces are applied only to a |
| 155 |
> |
pre-determined number of the outermost molecules. |
| 156 |
|
|
| 157 |
+ |
Beglov and Roux have developed a boundary model in which the hard |
| 158 |
+ |
sphere boundary has a radius that varies with the instantaneous |
| 159 |
+ |
configuration of the solute (and solvent) molecules.\cite{beglov:9050} |
| 160 |
+ |
This model contains a clear pressure and surface tension contribution |
| 161 |
+ |
to the free energy which XXX. |
| 162 |
+ |
|
| 163 |
+ |
\subsection*{Restraining Potentials} |
| 164 |
|
Restraining {\it potentials} introduce repulsive potentials at the |
| 165 |
|
surface of a sphere or other geometry. The solute and any explicit |
| 166 |
< |
solvent are therefore restrained inside this potential. Often the |
| 167 |
< |
potentials include a weak short-range attraction to maintain the |
| 168 |
< |
correct density at the boundary. Beglov and Roux have also introduced |
| 169 |
< |
a restraining boundary potential which relaxes dynamically depending |
| 170 |
< |
on the solute geometry and the force the explicit system exerts on the |
| 171 |
< |
shell.\cite{Beglov:1995fk} |
| 166 |
> |
solvent are therefore restrained inside the range defined by the |
| 167 |
> |
external potential. Often the potentials include a weak short-range |
| 168 |
> |
attraction to maintain the correct density at the boundary. Beglov |
| 169 |
> |
and Roux have also introduced a restraining boundary potential which |
| 170 |
> |
relaxes dynamically depending on the solute geometry and the force the |
| 171 |
> |
explicit system exerts on the shell.\cite{Beglov:1995fk} |
| 172 |
|
|
| 173 |
< |
Recently, Krilov {\it et al.} introduced a flexible boundary model |
| 174 |
< |
that uses a Lennard-Jones potential between the solvent molecules and |
| 175 |
< |
a boundary which is determined dynamically from the position of the |
| 176 |
< |
nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This approach allows |
| 177 |
< |
the confining potential to prevent solvent molecules from migrating |
| 178 |
< |
too far from the solute surface, while providing a weak attractive |
| 179 |
< |
force pulling the solvent molecules towards a fictitious bulk solvent. |
| 180 |
< |
Although this approach is appealing and has physical motivation, |
| 181 |
< |
nanoparticles do not deform far from their original geometries even at |
| 182 |
< |
temperatures which vaporize the nearby solvent. For the systems like |
| 183 |
< |
the one described, the flexible boundary model will be nearly |
| 173 |
> |
Recently, Krilov {\it et al.} introduced a {\it flexible} boundary |
| 174 |
> |
model that uses a Lennard-Jones potential between the solvent |
| 175 |
> |
molecules and a boundary which is determined dynamically from the |
| 176 |
> |
position of the nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This |
| 177 |
> |
approach allows the confining potential to prevent solvent molecules |
| 178 |
> |
from migrating too far from the solute surface, while providing a weak |
| 179 |
> |
attractive force pulling the solvent molecules towards a fictitious |
| 180 |
> |
bulk solvent. Although this approach is appealing and has physical |
| 181 |
> |
motivation, nanoparticles do not deform far from their original |
| 182 |
> |
geometries even at temperatures which vaporize the nearby solvent. For |
| 183 |
> |
the systems like this, the flexible boundary model will be nearly |
| 184 |
|
identical to a fixed-volume restraining potential. |
| 185 |
|
|
| 186 |
+ |
\subsection*{Hull methods} |
| 187 |
|
The approach of Kohanoff, Caro, and Finnis is the most promising of |
| 188 |
|
the methods for introducing both constant pressure and temperature |
| 189 |
|
into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru} |
| 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. Section \ref{sec:meth} |
| 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. |
| 212 |
|
|
| 213 |
|
\section{Methodology} |
| 214 |
|
\label{sec:meth} |
| 215 |
|
|
| 216 |
< |
We have developed a new method which uses an external bath at a fixed |
| 217 |
< |
constant pressure ($P$) and temperature ($T$). This bath interacts |
| 218 |
< |
only with the objects on the exterior hull of the system. Defining |
| 219 |
< |
the hull of the simulation is done in a manner similar to the approach |
| 220 |
< |
of Kohanoff, Caro and Finnis.\cite{Kohanoff:2005qm} That is, any |
| 221 |
< |
instantaneous configuration of the atoms in the system is considered |
| 222 |
< |
as a point cloud in three dimensional space. Delaunay triangulation |
| 223 |
< |
is used to find all facets between coplanar neighbors.\cite{DELAUNAY} |
| 224 |
< |
In highly symmetric point clouds, facets can contain many atoms, but |
| 225 |
< |
in all but the most symmetric of cases the facets are simple triangles |
| 226 |
< |
in 3-space that contain exactly three atoms. |
| 216 |
> |
The Langevin Hull uses an external bath at a fixed constant pressure |
| 217 |
> |
($P$) and temperature ($T$). This bath interacts only with the |
| 218 |
> |
objects on the exterior hull of the system. Defining the hull of the |
| 219 |
> |
simulation is done in a manner similar to the approach of Kohanoff, |
| 220 |
> |
Caro and Finnis.\cite{Kohanoff:2005qm} That is, any instantaneous |
| 221 |
> |
configuration of the atoms in the system is considered as a point |
| 222 |
> |
cloud in three dimensional space. Delaunay triangulation is used to |
| 223 |
> |
find all facets between coplanar |
| 224 |
> |
neighbors.\cite{delaunay,springerlink:10.1007/BF00977785} In highly |
| 225 |
> |
symmetric point clouds, facets can contain many atoms, but in all but |
| 226 |
> |
the most symmetric of cases the facets are simple triangles in 3-space |
| 227 |
> |
that contain exactly three atoms. |
| 228 |
|
|
| 229 |
|
The convex hull is the set of facets that have {\it no concave |
| 230 |
< |
corners} at an atomic site. This eliminates all facets on the |
| 231 |
< |
interior of the point cloud, leaving only those exposed to the |
| 232 |
< |
bath. Sites on the convex hull are dynamic. As molecules re-enter the |
| 233 |
< |
cluster, all interactions between atoms on that molecule and the |
| 234 |
< |
external bath are removed. Since the edge is determined dynamically |
| 235 |
< |
as the simulation progresses, no {\it a priori} geometry is |
| 236 |
< |
defined. The pressure and temperature bath interacts {\it directly} |
| 230 |
> |
corners} at an atomic site.\cite{Barber96,EDELSBRUNNER:1994oq} This |
| 231 |
> |
eliminates all facets on the interior of the point cloud, leaving only |
| 232 |
> |
those exposed to the bath. Sites on the convex hull are dynamic; as |
| 233 |
> |
molecules re-enter the cluster, all interactions between atoms on that |
| 234 |
> |
molecule and the external bath are removed. Since the edge is |
| 235 |
> |
determined dynamically as the simulation progresses, no {\it a priori} |
| 236 |
> |
geometry is defined. The pressure and temperature bath interacts only |
| 237 |
|
with the atoms on the edge and not with atoms interior to the |
| 238 |
|
simulation. |
| 239 |
|
|
| 221 |
– |
|
| 240 |
|
\begin{figure} |
| 241 |
|
\includegraphics[width=\linewidth]{hullSample} |
| 242 |
|
\caption{The external temperature and pressure bath interacts only |
| 246 |
|
\label{fig:hullSample} |
| 247 |
|
\end{figure} |
| 248 |
|
|
| 249 |
< |
|
| 232 |
< |
Atomic sites in the interior of the point cloud move under standard |
| 249 |
> |
Atomic sites in the interior of the simulation move under standard |
| 250 |
|
Newtonian dynamics, |
| 251 |
|
\begin{equation} |
| 252 |
|
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
| 262 |
|
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
| 263 |
|
\end{equation} |
| 264 |
|
|
| 265 |
< |
The external bath interacts directly with the facets of the convex |
| 266 |
< |
hull. Since each vertex (or atom) provides one corner of a triangular |
| 267 |
< |
facet, the force on the facets are divided equally to each vertex. |
| 268 |
< |
However, each vertex can participate in multiple facets, so the resultant |
| 269 |
< |
force is a sum over all facets $f$ containing vertex $i$: |
| 265 |
> |
The external bath interacts indirectly with the atomic sites through |
| 266 |
> |
the intermediary of the hull facets. Since each vertex (or atom) |
| 267 |
> |
provides one corner of a triangular facet, the force on the facets are |
| 268 |
> |
divided equally to each vertex. However, each vertex can participate |
| 269 |
> |
in multiple facets, so the resultant force is a sum over all facets |
| 270 |
> |
$f$ containing vertex $i$: |
| 271 |
|
\begin{equation} |
| 272 |
|
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\ |
| 273 |
|
} f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf |
| 285 |
|
& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t) |
| 286 |
|
\end{array} |
| 287 |
|
\end{equation} |
| 288 |
< |
Here, $A_f$ and $\hat{n}_f$ are the area and normal vectors for facet |
| 289 |
< |
$f$, respectively. ${\mathbf v}_f(t)$ is the velocity of the facet |
| 290 |
< |
centroid, |
| 288 |
> |
Here, $A_f$ and $\hat{n}_f$ are the area and (outward-facing) normal |
| 289 |
> |
vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is the |
| 290 |
> |
velocity of the facet centroid, |
| 291 |
|
\begin{equation} |
| 292 |
|
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
| 293 |
|
\end{equation} |
| 303 |
|
\label{eq:randomForce} |
| 304 |
|
\end{eqnarray} |
| 305 |
|
|
| 306 |
< |
Once the resistance tensor is known for a given facet a stochastic |
| 306 |
> |
Once the resistance tensor is known for a given facet, a stochastic |
| 307 |
|
vector that has the properties in Eq. (\ref{eq:randomForce}) can be |
| 308 |
< |
done efficiently by carrying out a Cholesky decomposition to obtain |
| 309 |
< |
the square root matrix of the resistance tensor, |
| 308 |
> |
calculated efficiently by carrying out a Cholesky decomposition to |
| 309 |
> |
obtain the square root matrix of the resistance tensor, |
| 310 |
|
\begin{equation} |
| 311 |
|
\Xi_f = {\bf S} {\bf S}^{T}, |
| 312 |
|
\label{eq:Cholesky} |
| 333 |
|
prohibitively expensive if it were recomputed at each step in a |
| 334 |
|
molecular dynamics simulation. |
| 335 |
|
|
| 336 |
< |
We are utilizing an approximate resistance tensor obtained by first |
| 337 |
< |
constructing the Oseen tensor for the interaction of the centroid of |
| 338 |
< |
the facet ($f$) with each of the subfacets $j$, |
| 336 |
> |
Instead, we are utilizing an approximate resistance tensor obtained by |
| 337 |
> |
first constructing the Oseen tensor for the interaction of the |
| 338 |
> |
centroid of the facet ($f$) with each of the subfacets $\ell=1,2,3$, |
| 339 |
|
\begin{equation} |
| 340 |
< |
T_{jf}=\frac{A_j}{8\pi\eta R_{jf}}\left(I + |
| 341 |
< |
\frac{\mathbf{R}_{jf}\mathbf{R}_{jf}^T}{R_{jf}^2}\right) |
| 340 |
> |
T_{\ell f}=\frac{A_\ell}{8\pi\eta R_{\ell f}}\left(I + |
| 341 |
> |
\frac{\mathbf{R}_{\ell f}\mathbf{R}_{\ell f}^T}{R_{\ell f}^2}\right) |
| 342 |
|
\end{equation} |
| 343 |
< |
Here, $A_j$ is the area of subfacet $j$ which is a triangle containing |
| 344 |
< |
two of the vertices of the facet along with the centroid. |
| 345 |
< |
$\mathbf{R}_{jf}$ is the vector between the centroid of facet $f$ and |
| 346 |
< |
the centroid of sub-facet $j$, and $I$ is the ($3 \times 3$) identity |
| 347 |
< |
matrix. $\eta$ is the viscosity of the external bath. |
| 343 |
> |
Here, $A_\ell$ is the area of subfacet $\ell$ which is a triangle |
| 344 |
> |
containing two of the vertices of the facet along with the centroid. |
| 345 |
> |
$\mathbf{R}_{\ell f}$ is the vector between the centroid of facet $f$ |
| 346 |
> |
and the centroid of sub-facet $\ell$, and $I$ is the ($3 \times 3$) |
| 347 |
> |
identity matrix. $\eta$ is the viscosity of the external bath. |
| 348 |
|
|
| 349 |
|
\begin{figure} |
| 350 |
|
\includegraphics[width=\linewidth]{hydro} |
| 357 |
|
\label{hydro} |
| 358 |
|
\end{figure} |
| 359 |
|
|
| 360 |
< |
The Oseen tensors for each of the sub-facets are added together, and |
| 361 |
< |
the resulting matrix is inverted to give a $3 \times 3$ resistance |
| 362 |
< |
tensor for translations of the triangular facet, |
| 360 |
> |
The tensors for each of the sub-facets are added together, and the |
| 361 |
> |
resulting matrix is inverted to give a $3 \times 3$ resistance tensor |
| 362 |
> |
for translations of the triangular facet, |
| 363 |
|
\begin{equation} |
| 364 |
|
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
| 365 |
|
\end{equation} |
| 369 |
|
configurations, so this appears to be a reasonably good approximation. |
| 370 |
|
|
| 371 |
|
We have implemented this method by extending the Langevin dynamics |
| 372 |
< |
integrator in our code, OpenMD.\cite{Meineke2005,openmd} The Delaunay |
| 373 |
< |
triangulation and computation of the convex hull are done using calls |
| 374 |
< |
to the qhull library.\cite{Qhull} There is a moderate penalty for |
| 375 |
< |
computing the convex hull at each step in the molecular dynamics |
| 376 |
< |
simulation (HOW MUCH?), but the convex hull is remarkably easy to |
| 377 |
< |
parallelize on distributed memory machines (see Appendix A). |
| 372 |
> |
integrator in our code, OpenMD.\cite{Meineke2005,openmd} At each |
| 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 |
| 377 |
> |
configuration. |
| 378 |
> |
\item The convex hull is computed and facets are identified. |
| 379 |
> |
\item For each facet: |
| 380 |
> |
\begin{itemize} |
| 381 |
> |
\item[a.] The force from the pressure bath ($-PA_f\hat{n}_f$) is |
| 382 |
> |
computed. |
| 383 |
> |
\item[b.] The resistance tensor ($\Xi_f(t)$) is computed using the |
| 384 |
> |
viscosity ($\eta$) of the bath. |
| 385 |
> |
\item[c.] Facet drag ($-\Xi_f(t) \mathbf{v}_f(t)$) forces are |
| 386 |
> |
computed. |
| 387 |
> |
\item[d.] Random forces ($\mathbf{R}_f(t)$) are computed using the |
| 388 |
> |
resistance tensor and the temperature ($T$) of the bath. |
| 389 |
> |
\end{itemize} |
| 390 |
> |
\item The facet forces are divided equally among the vertex atoms. |
| 391 |
> |
\item Atomic positions and velocities are propagated. |
| 392 |
> |
\end{enumerate} |
| 393 |
> |
The Delaunay triangulation and computation of the convex hull are done |
| 394 |
> |
using calls to the qhull library.\cite{Qhull} There is a minimal |
| 395 |
> |
penalty for computing the convex hull and resistance tensors at each |
| 396 |
> |
step in the molecular dynamics simulation (roughly 0.02 $\times$ cost |
| 397 |
> |
of a single force evaluation), and the convex hull is remarkably easy |
| 398 |
> |
to parallelize on distributed memory machines (see Appendix A). |
| 399 |
|
|
| 400 |
|
\section{Tests \& Applications} |
| 401 |
|
\label{sec:tests} |
| 402 |
|
|
| 403 |
|
To test the new method, we have carried out simulations using the |
| 404 |
|
Langevin Hull on: 1) a crystalline system (gold nanoparticles), 2) a |
| 405 |
< |
liquid droplet (SPC/E water),\cite{SPCE} and 3) a heterogeneous |
| 406 |
< |
mixture (gold nanoparticles in a water droplet). In each case, we have |
| 407 |
< |
computed properties that depend on the external applied pressure. Of |
| 408 |
< |
particular interest for the single-phase systems is the bulk modulus, |
| 405 |
> |
liquid droplet (SPC/E water),\cite{Berendsen1987} and 3) a |
| 406 |
> |
heterogeneous mixture (gold nanoparticles in a water droplet). In each |
| 407 |
> |
case, we have computed properties that depend on the external applied |
| 408 |
> |
pressure. Of particular interest for the single-phase systems is the |
| 409 |
> |
isothermal compressibility, |
| 410 |
|
\begin{equation} |
| 411 |
|
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right |
| 412 |
|
)_{T}. |
| 417 |
|
simulation boxes is that the volume of a three-dimensional point cloud |
| 418 |
|
is not well-defined. In order to compute the compressibility of a |
| 419 |
|
bulk material, we make an assumption that the number density, $\rho = |
| 420 |
< |
\frac{N}{V}$, is uniform within some region of the cloud. The |
| 420 |
> |
\frac{N}{V}$, is uniform within some region of the point cloud. The |
| 421 |
|
compressibility can then be expressed in terms of the average number |
| 422 |
|
of particles in that region, |
| 423 |
|
\begin{equation} |
| 424 |
< |
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right |
| 424 |
> |
\kappa_{T} = -\frac{1}{N} \left ( \frac{\partial N}{\partial P} \right |
| 425 |
|
)_{T} |
| 426 |
|
\label{eq:BMN} |
| 427 |
|
\end{equation} |
| 429 |
|
the middle of the cluster. $N$ is the average number of molecules |
| 430 |
|
found within this region throughout a given simulation. The geometry |
| 431 |
|
and size of the region is arbitrary, and any bulk-like portion of the |
| 432 |
< |
cluster can be used to compute the bulk modulus. |
| 432 |
> |
cluster can be used to compute the compressibility. |
| 433 |
|
|
| 434 |
< |
One might assume that the volume of the convex hull could be taken as |
| 435 |
< |
the system volume in the compressibility expression (Eq. \ref{eq:BM}), |
| 436 |
< |
but this has implications at lower pressures (which are explored in |
| 437 |
< |
detail in the section on water droplets). |
| 434 |
> |
One might assume that the volume of the convex hull could simply be |
| 435 |
> |
taken as the system volume $V$ in the compressibility expression |
| 436 |
> |
(Eq. \ref{eq:BM}), but this has implications at lower pressures (which |
| 437 |
> |
are explored in detail in the section on water droplets). |
| 438 |
|
|
| 439 |
|
The metallic force field in use for the gold nanoparticles is the |
| 440 |
|
quantum Sutton-Chen (QSC) model.\cite{PhysRevB.59.3527} In all |
| 449 |
|
|
| 450 |
|
\subsection{Bulk modulus of gold nanoparticles} |
| 451 |
|
|
| 452 |
< |
The bulk modulus is well-known for gold, and it provides a good first |
| 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. |
| 454 |
|
|
| 415 |
– |
|
| 455 |
|
\begin{figure} |
| 456 |
< |
\includegraphics[width=\linewidth]{pressure_tb} |
| 457 |
< |
\caption{Pressure response is rapid (18 \AA gold nanoparticle), target |
| 458 |
< |
pressure = 4 GPa} |
| 456 |
> |
\includegraphics[width=\linewidth]{P_T_combined} |
| 457 |
> |
\caption{Pressure and temperature response of an 18 \AA\ gold |
| 458 |
> |
nanoparticle initially when first placed in the Langevin Hull |
| 459 |
> |
($T_\mathrm{bath}$ = 300K, $P_\mathrm{bath}$ = 4 GPa) and starting |
| 460 |
> |
from initial conditions that were far from the bath pressure and |
| 461 |
> |
temperature. The pressure response is rapid, and the thermal |
| 462 |
> |
equilibration depends on both total surface area and the viscosity |
| 463 |
> |
of the bath.} |
| 464 |
|
\label{pressureResponse} |
| 421 |
– |
\end{figure} |
| 422 |
– |
|
| 423 |
– |
\begin{figure} |
| 424 |
– |
\includegraphics[width=\linewidth]{temperature_tb} |
| 425 |
– |
\caption{Temperature equilibration depends on surface area and bath |
| 426 |
– |
viscosity. Target Temperature = 300K} |
| 427 |
– |
\label{temperatureResponse} |
| 465 |
|
\end{figure} |
| 466 |
|
|
| 467 |
|
\begin{equation} |
| 511 |
|
\includegraphics[width=\linewidth]{flytest2} |
| 512 |
|
\caption{At low pressures, the liquid is in equilibrium with the vapor |
| 513 |
|
phase, and isolated molecules can detach from the liquid droplet. |
| 514 |
< |
This is expected behavior, but the reported volume of the convex |
| 515 |
< |
hull includes large regions of empty space. For this reason, |
| 514 |
> |
This is expected behavior, but the volume of the convex hull |
| 515 |
> |
includes large regions of empty space. For this reason, |
| 516 |
|
compressibilities are computed using local number densities rather |
| 517 |
|
than hull volumes.} |
| 518 |
|
\label{fig:coneOfShame} |
| 519 |
|
\end{figure} |
| 520 |
|
|
| 521 |
< |
At higher pressures, the equilibrium favors the liquid phase, and the |
| 522 |
< |
hull geometries are much more compact. Because of the liquid-vapor |
| 523 |
< |
effect on the convex hull, the regional number density approach |
| 524 |
< |
(Eq. \ref{eq:BMN}) provides more reliable estimates of the bulk |
| 525 |
< |
modulus. |
| 521 |
> |
At higher pressures, the equilibrium strongly favors the liquid phase, |
| 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. |
| 526 |
|
|
| 527 |
< |
We initially used the classic compressibility formula to calculate the the isothermal compressibility at each target pressure. These calculations yielded compressibility values that were dramatically higher than both previous simulations and experiment. The particular compressibility expression used requires the calculation of both a volume and pressure differential, thereby stipulating that the data from at least two simulations at different pressures must be used to calculate the isothermal compressibility at one pressure. |
| 528 |
< |
|
| 529 |
< |
Regardless of the difficulty in obtaining accurate hull |
| 530 |
< |
volumes at low temperature and pressures, the Langevin Hull NPT method |
| 531 |
< |
provides reasonable isothermal compressibility values for water |
| 532 |
< |
through a large range of pressures. |
| 496 |
< |
|
| 497 |
< |
Per the fluctuation dissipation theorem \cite{Debenedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure |
| 498 |
< |
|
| 527 |
> |
In both the traditional compressibility formula (Eq. \ref{eq:BM}) and |
| 528 |
> |
the number density version (Eq. \ref{eq:BMN}), multiple simulations at |
| 529 |
> |
different pressures must be done to compute the first derivatives. It |
| 530 |
> |
is also possible to compute the compressibility using the fluctuation |
| 531 |
> |
dissipation theorem using either fluctuations in the |
| 532 |
> |
volume,\cite{Debenedetti1986}, |
| 533 |
|
\begin{equation} |
| 534 |
< |
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T} |
| 534 |
> |
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle |
| 535 |
> |
V \right \rangle ^{2}}{V \, k_{B} \, T}, |
| 536 |
|
\end{equation} |
| 537 |
+ |
or, equivalently, fluctuations in the number of molecules within the |
| 538 |
+ |
fixed region, |
| 539 |
+ |
\begin{equation} |
| 540 |
+ |
\kappa_{T} = \frac{\left \langle N^{2} \right \rangle - \left \langle |
| 541 |
+ |
N \right \rangle ^{2}}{N \, k_{B} \, T}, |
| 542 |
+ |
\end{equation} |
| 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 |
| 547 |
+ |
calculation that relies on the hull volume will suffer these effects. |
| 548 |
+ |
WE NEED MORE HERE. |
| 549 |
|
|
| 503 |
– |
Thus, the compressibility of each simulation run can be calculated entirely independently from all other trajectories. However, the resulting compressibilities were still as much as an order of magnitude larger than the reference values. The effect was particularly pronounced at the low end of the pressure range. At ambient temperature and low pressures, there exists an equilibrium between vapor and liquid phases. Vapor molecules are naturally more diffuse around the exterior of the cluster, causing artificially large cluster volumes. Any compressibility calculation that relies on the hull volume will suffer these effects. |
| 504 |
– |
|
| 505 |
– |
|
| 550 |
|
\subsection{Molecular orientation distribution at cluster boundary} |
| 551 |
|
|
| 552 |
< |
In order for non-periodic boundary conditions to be widely applicable, they must be constructed in such a way that they allow a finite, usually small, simulated system to replicate the properties of an infinite bulk system. Naturally, this requirement has spawned many methods for inserting boundaries into simulated systems [REF... ?]. Of particular interest to our characterization of the Langevin Hull is the orientation of water molecules included in the geometric hull. Ideally, all molecules in the cluster will have the same orientational distribution as bulk water. |
| 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. |
| 560 |
|
|
| 561 |
< |
The orientation of molecules at the edges of a simulated cluster has long been a concern when performing simulations of explicitly non-periodic systems. Early work led to the surface constrained soft sphere dipole model (SCSSD) \cite{Warshel1978} in which the surface molecules are fixed in a random orientation representative of the bulk solvent structural properties. Belch, et al \cite{Belch1985} simulated clusters of TIPS2 water surrounded by a hydrophobic bounding potential. The spherical hydrophobic boundary induced dangling hydrogen bonds at the surface that propagated deep into the cluster, affecting 70\% of the 100 molecules in the simulation. This result echoes an earlier study which showed that an extended planar hydrophobic surface caused orientational preference at the surface which extended 7 \r{A} into the liquid simulation cell \cite{Lee1984}. The surface constrained all-atom solvent (SCAAS) model \cite{King1989} improved upon its SCSSD predecessor. The SCAAS model utilizes a polarization constraint which is applied to the surface molecules to maintain bulk-like structure at the cluster surface. A radial constraint is used to maintain the desired bulk density of the liquid. Both constraint forces are applied only to a pre-determined number of the outermost molecules. |
| 561 |
> |
The orientation of molecules at the edges of a simulated cluster has |
| 562 |
> |
long been a concern when performing simulations of explicitly |
| 563 |
> |
non-periodic systems. Early work led to the surface constrained soft |
| 564 |
> |
sphere dipole model (SCSSD) \cite{Warshel1978} in which the surface |
| 565 |
> |
molecules are fixed in a random orientation representative of the bulk |
| 566 |
> |
solvent structural properties. Belch, et al \cite{Belch1985} simulated |
| 567 |
> |
clusters of TIPS2 water surrounded by a hydrophobic bounding |
| 568 |
> |
potential. The spherical hydrophobic boundary induced dangling |
| 569 |
> |
hydrogen bonds at the surface that propagated deep into the cluster, |
| 570 |
> |
affecting 70\% of the 100 molecules in the simulation. This result |
| 571 |
> |
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 of molecules be fixed, nor does it utilize an explicitly hydrophobic boundary, orientational constraint or radial constraint. The number and identity of the molecules included on the convex hull are dynamic properties, thus avoiding the formation of an artificial solvent boundary layer. The hope is that the water molecules on the surface of the cluster, if left to their own devices in the absence of orientational and radial constraints, will maintain a bulk-like orientational distribution. |
| 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. |
| 593 |
|
|
| 621 |
|
|
| 622 |
|
\subsection{Heterogeneous nanoparticle / water mixtures} |
| 623 |
|
|
| 624 |
+ |
\section{Discussion} |
| 625 |
+ |
\label{sec:discussion} |
| 626 |
|
|
| 627 |
|
\section*{Appendix A: Computing Convex Hulls on Parallel Computers} |
| 628 |
|
|
