| 195 |
|
\section{Methodology} |
| 196 |
|
\label{sec:meth} |
| 197 |
|
|
| 198 |
< |
We have developed a new method which uses a constant pressure and |
| 199 |
< |
temperature bath. This bath interacts only with the objects that are |
| 200 |
< |
currently at the edge of the system. Since the edge is determined |
| 201 |
< |
dynamically as the simulation progresses, no {\it a priori} geometry |
| 202 |
< |
is defined. The pressure and temperature bath interacts {\it |
| 203 |
< |
directly} with the atoms on the edge and not with atoms interior to |
| 204 |
< |
the simulation. This means that there are no affine transforms |
| 205 |
< |
required. There are also no fictitious particles or bounding |
| 206 |
< |
potentials used in this approach. |
| 207 |
< |
|
| 208 |
< |
The basics of the method are as follows. The simulation starts as a |
| 209 |
< |
collection of atomic locations in three dimensions (a point cloud). |
| 210 |
< |
Delaunay triangulation is used to find all facets between coplanar |
| 211 |
< |
neighbors. In highly symmetric point clouds, facets can contain many |
| 212 |
< |
atoms, but in all but the most symmetric of cases one might experience |
| 213 |
< |
in a molecular dynamics simulation, the facets are simple triangles in |
| 214 |
< |
3-space that contain exactly three atoms. |
| 198 |
> |
We have developed a new method which uses an external bath at a fixed |
| 199 |
> |
constant pressure ($P$) and temperature ($T$). This bath interacts |
| 200 |
> |
only with the objects on the exterior hull of the system. Defining |
| 201 |
> |
the hull of the simulation is done in a manner similar to the approach |
| 202 |
> |
of Kohanoff, Caro and Finnis.\cite{Kohanoff:2005qm} That is, any |
| 203 |
> |
instantaneous configuration of the atoms in the system is considered |
| 204 |
> |
as a point cloud in three dimensional space. Delaunay triangulation |
| 205 |
> |
is used to find all facets between coplanar neighbors.\cite{DELAUNAY} |
| 206 |
> |
In highly symmetric point clouds, facets can contain many atoms, but |
| 207 |
> |
in all but the most symmetric of cases the facets are simple triangles |
| 208 |
> |
in 3-space that contain exactly three atoms. |
| 209 |
|
|
| 210 |
|
The convex hull is the set of facets that have {\it no concave |
| 211 |
|
corners} at an atomic site. This eliminates all facets on the |
| 212 |
|
interior of the point cloud, leaving only those exposed to the |
| 213 |
|
bath. Sites on the convex hull are dynamic. As molecules re-enter the |
| 214 |
|
cluster, all interactions between atoms on that molecule and the |
| 215 |
< |
external bath are removed. |
| 215 |
> |
external bath are removed. Since the edge is determined dynamically |
| 216 |
> |
as the simulation progresses, no {\it a priori} geometry is |
| 217 |
> |
defined. The pressure and temperature bath interacts {\it directly} |
| 218 |
> |
with the atoms on the edge and not with atoms interior to the |
| 219 |
> |
simulation. |
| 220 |
|
|
| 221 |
< |
For atomic sites in the interior of the point cloud, the equations of |
| 222 |
< |
motion are simple Newtonian dynamics, |
| 221 |
> |
Atomic sites in the interior of the point cloud move under standard |
| 222 |
> |
Newtonian dynamics, |
| 223 |
|
\begin{equation} |
| 224 |
|
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
| 225 |
|
\label{eq:Newton} |
| 235 |
|
\end{equation} |
| 236 |
|
|
| 237 |
|
The external bath interacts directly with the facets of the convex |
| 238 |
< |
hull. Since each vertex (or atom) provides one corner of a triangular |
| 238 |
> |
hull. Since each vertex (or atom) provides one corner of a triangular |
| 239 |
|
facet, the force on the facets are divided equally to each vertex. |
| 240 |
|
However, each vertex can participate in multiple facets, so the resultant |
| 241 |
|
force is a sum over all facets $f$ containing vertex $i$: |
| 247 |
|
|
| 248 |
|
The external pressure bath applies a force to the facets of the convex |
| 249 |
|
hull in direct proportion to the area of the facet, while the thermal |
| 250 |
< |
coupling depends on the solvent temperature, friction and the size and |
| 251 |
< |
shape of each facet. The thermal interactions are expressed as a |
| 252 |
< |
typical Langevin description of the forces, |
| 250 |
> |
coupling depends on the solvent temperature, viscosity and the size |
| 251 |
> |
and shape of each facet. The thermal interactions are expressed as a |
| 252 |
> |
standard Langevin description of the forces, |
| 253 |
|
\begin{equation} |
| 254 |
|
\begin{array}{rclclcl} |
| 255 |
|
{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\ |
| 256 |
|
& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t) |
| 257 |
|
\end{array} |
| 258 |
|
\end{equation} |
| 259 |
< |
Here, $P$ is the external pressure, $A_f$ and $\hat{n}_f$ are the area |
| 260 |
< |
and normal vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is |
| 261 |
< |
the velocity of the facet, |
| 259 |
> |
Here, $A_f$ and $\hat{n}_f$ are the area and normal vectors for facet |
| 260 |
> |
$f$, respectively. ${\mathbf v}_f(t)$ is the velocity of the facet |
| 261 |
> |
centroid, |
| 262 |
|
\begin{equation} |
| 263 |
|
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
| 264 |
|
\end{equation} |
| 265 |
< |
and $\Xi_f(t)$ is an approximate ($3 \times 3$) hydrodynamic tensor |
| 266 |
< |
that depends on the geometry and surface area of facet $f$ and the |
| 267 |
< |
viscosity of the fluid (See Appendix A). The hydrodynamic tensor is |
| 268 |
< |
related to the fluctuations of the random force, $\mathbf{R}(t)$, by |
| 269 |
< |
the fluctuation-dissipation theorem, |
| 265 |
> |
and $\Xi_f(t)$ is an approximate ($3 \times 3$) resistance tensor that |
| 266 |
> |
depends on the geometry and surface area of facet $f$ and the |
| 267 |
> |
viscosity of the fluid. The resistance tensor is related to the |
| 268 |
> |
fluctuations of the random force, $\mathbf{R}(t)$, by the |
| 269 |
> |
fluctuation-dissipation theorem, |
| 270 |
|
\begin{eqnarray} |
| 271 |
|
\left< {\mathbf R}_f(t) \right> & = & 0 \\ |
| 272 |
|
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\ |
| 274 |
|
\label{eq:randomForce} |
| 275 |
|
\end{eqnarray} |
| 276 |
|
|
| 277 |
< |
Once the hydrodynamic tensor is known for a given facet (see Appendix |
| 278 |
< |
A) obtaining a stochastic vector that has the properties in |
| 279 |
< |
Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a |
| 280 |
< |
one-time Cholesky decomposition to obtain the square root matrix of |
| 283 |
< |
the resistance tensor, |
| 277 |
> |
Once the resistance tensor is known for a given facet a stochastic |
| 278 |
> |
vector that has the properties in Eq. (\ref{eq:randomForce}) can be |
| 279 |
> |
done efficiently by carrying out a Cholesky decomposition to obtain |
| 280 |
> |
the square root matrix of the resistance tensor, |
| 281 |
|
\begin{equation} |
| 282 |
|
\Xi_f = {\bf S} {\bf S}^{T}, |
| 283 |
|
\label{eq:Cholesky} |
| 294 |
|
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
| 295 |
|
have the correct properties required by Eq. (\ref{eq:randomForce}). |
| 296 |
|
|
| 297 |
< |
Our treatment of the hydrodynamic tensor must be approximate. $\Xi$ |
| 298 |
< |
for a triangular plate would normally be treated as a $6 \times 6$ |
| 297 |
> |
Our treatment of the resistance tensor is approximate. $\Xi$ for a |
| 298 |
> |
rigid triangular plate would normally be treated as a $6 \times 6$ |
| 299 |
|
tensor that includes translational and rotational drag as well as |
| 300 |
< |
translational-rotational coupling. The computation of hydrodynamic |
| 300 |
> |
translational-rotational coupling. The computation of resistance |
| 301 |
|
tensors for rigid bodies has been detailed |
| 302 |
|
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun2008} |
| 303 |
|
but the standard approach involving bead approximations would be |
| 304 |
|
prohibitively expensive if it were recomputed at each step in a |
| 305 |
|
molecular dynamics simulation. |
| 306 |
|
|
| 307 |
< |
We are utilizing an approximate hydrodynamic tensor obtained by first |
| 307 |
> |
We are utilizing an approximate resistance tensor obtained by first |
| 308 |
|
constructing the Oseen tensor for the interaction of the centroid of |
| 309 |
|
the facet ($f$) with each of the subfacets $j$, |
| 310 |
|
\begin{equation} |
| 319 |
|
|
| 320 |
|
\begin{figure} |
| 321 |
|
\includegraphics[width=\linewidth]{hydro} |
| 322 |
< |
\caption{The hydrodynamic tensor $\Xi$ for a facet comprising sites $i$, |
| 323 |
< |
$j$, and $k$ is constructed using Oseen tensor contributions |
| 324 |
< |
between the centoid of the facet $f$ and each of the sub-facets |
| 325 |
< |
($i,f,j$), ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets |
| 326 |
< |
are located at $1$, $2$, and $3$, and the area of each sub-facet is |
| 322 |
> |
\caption{The resistance tensor $\Xi$ for a facet comprising sites $i$, |
| 323 |
> |
$j$, and $k$ is constructed using Oseen tensor contributions between |
| 324 |
> |
the centoid of the facet $f$ and each of the sub-facets ($i,f,j$), |
| 325 |
> |
($j,f,k$), and ($k,f,i$). The centroids of the sub-facets are |
| 326 |
> |
located at $1$, $2$, and $3$, and the area of each sub-facet is |
| 327 |
|
easily computed using half the cross product of two of the edges.} |
| 328 |
|
\label{hydro} |
| 329 |
|
\end{figure} |
| 330 |
|
|
| 331 |
< |
The Oseen tensors for each of the sub-facets are summed, and the |
| 332 |
< |
resulting matrix is inverted to give a $3 \times 3$ hydrodynamic |
| 333 |
< |
tensor for translations of the triangular plate, |
| 331 |
> |
The Oseen tensors for each of the sub-facets are added together, and |
| 332 |
> |
the resulting matrix is inverted to give a $3 \times 3$ resistance |
| 333 |
> |
tensor for translations of the triangular facet, |
| 334 |
|
\begin{equation} |
| 335 |
|
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
| 336 |
|
\end{equation} |
| 337 |
+ |
Note that this treatment explicitly ignores rotations (and |
| 338 |
+ |
translational-rotational coupling) of the facet. In compact systems, |
| 339 |
+ |
the facets stay relatively fixed in orientation between |
| 340 |
+ |
configurations, so this appears to be a reasonably good approximation. |
| 341 |
+ |
|
| 342 |
|
We have implemented this method by extending the Langevin dynamics |
| 343 |
< |
integrator in our group code, OpenMD.\cite{Meineke2005,openmd} There |
| 344 |
< |
is a moderate penalty for computing the convex hull at each step in |
| 345 |
< |
the molecular dynamics simulation (HOW MUCH?), but the convex hull is |
| 346 |
< |
remarkably easy to parallelize on distributed memory machines (see |
| 347 |
< |
Appendix B). |
| 343 |
> |
integrator in our code, OpenMD.\cite{Meineke2005,openmd} The Delaunay |
| 344 |
> |
triangulation and computation of the convex hull are done using calls |
| 345 |
> |
to the qhull library.\cite{qhull} There is a moderate penalty for |
| 346 |
> |
computing the convex hull at each step in the molecular dynamics |
| 347 |
> |
simulation (HOW MUCH?), but the convex hull is remarkably easy to |
| 348 |
> |
parallelize on distributed memory machines (see Appendix A). |
| 349 |
|
|
| 350 |
|
\section{Tests \& Applications} |
| 351 |
|
\label{sec:tests} |
| 352 |
|
|
| 353 |
+ |
In order to test this method, we have carried out simulations using |
| 354 |
+ |
the Langevin Hull on a crystalline system (gold nanoparticles), a |
| 355 |
+ |
liquid droplet (SPC/E water), and a heterogeneous mixture (gold |
| 356 |
+ |
nanoparticles in a water droplet). In each case, we have computed |
| 357 |
+ |
bulk properties that depend on the external applied pressure. Of |
| 358 |
+ |
particular interest is the bulk modulus, |
| 359 |
+ |
\begin{equation} |
| 360 |
+ |
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right |
| 361 |
+ |
)_{T}. |
| 362 |
+ |
\label{eq:BM} |
| 363 |
+ |
\end{equation} |
| 364 |
+ |
|
| 365 |
+ |
One problem with eliminating periodic boundary conditions and |
| 366 |
+ |
simulation boxes is that the volume of a three-dimensional point cloud |
| 367 |
+ |
is not well-defined. In order to compute the compressibility of a |
| 368 |
+ |
bulk material, we make an assumption that the number density, $\rho = |
| 369 |
+ |
\frac{N}{V}$, is uniform within some region of the cloud. The |
| 370 |
+ |
compressibility can then be expressed in terms of the average number |
| 371 |
+ |
of particles in that region, |
| 372 |
+ |
\begin{equation} |
| 373 |
+ |
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right |
| 374 |
+ |
)_{T} |
| 375 |
+ |
\label{eq:BMN} |
| 376 |
+ |
\end{equation} |
| 377 |
+ |
The region we pick is a spherical volume of 10 \AA radius centered in |
| 378 |
+ |
the middle of the cluster. This radius is arbitrary, and any |
| 379 |
+ |
bulk-like portion of the cluster can be used to compute the bulk |
| 380 |
+ |
modulus. |
| 381 |
+ |
|
| 382 |
+ |
One might also assume that the volume of the convex hull could be |
| 383 |
+ |
taken as the system volume in the compressibility expression |
| 384 |
+ |
(Eq. \ref{eq:BM}), but this has implications at lower pressures (which |
| 385 |
+ |
are explored in detail in the section on water droplets). |
| 386 |
+ |
|
| 387 |
|
\subsection{Bulk modulus of gold nanoparticles} |
| 388 |
|
|
| 389 |
|
\begin{figure} |
| 413 |
|
|
| 414 |
|
\subsection{Compressibility of SPC/E water clusters} |
| 415 |
|
|
| 416 |
< |
Both NVT \cite{Glattli2002} and NPT \cite{Motakabbir1990, Pi2009} molecular dynamics simulations of SPC/E water have yielded values for the isothermal compressibility of water that agree well with experiment \cite{Fine1973}. The results of three different methods for computing the isothermal compressibility from Langevin Hull simulations for pressures between 1 and 6500 atm are shown in Fig. 5 along with compressibility values obtained from both other SPC/E simulations and experiment. Compressibility values from all references are for applied pressures within the range 1 - 1000 atm. |
| 416 |
> |
Prior molecular dynamics simulations on SPC/E water (both in |
| 417 |
> |
NVT~\cite{Glattli2002} and NPT~\cite{Motakabbir1990, Pi2009} |
| 418 |
> |
ensembles) have yielded values for the isothermal compressibility that |
| 419 |
> |
agree well with experiment.\cite{Fine1973} The results of two |
| 420 |
> |
different approaches for computing the isothermal compressibility from |
| 421 |
> |
Langevin Hull simulations for pressures between 1 and 6500 atm are |
| 422 |
> |
shown in Fig. \ref{fig:compWater} along with compressibility values |
| 423 |
> |
obtained from both other SPC/E simulations and experiment. |
| 424 |
> |
Compressibility values from all references are for applied pressures |
| 425 |
> |
within the range 1 - 1000 atm. |
| 426 |
|
|
| 427 |
|
\begin{figure} |
| 428 |
< |
\includegraphics[width=\linewidth]{new_isothermal} |
| 428 |
> |
\includegraphics[width=\linewidth]{new_isothermalN} |
| 429 |
|
\caption{Compressibility of SPC/E water} |
| 430 |
< |
\label{compWater} |
| 430 |
> |
\label{fig:compWater} |
| 431 |
|
\end{figure} |
| 432 |
|
|
| 433 |
< |
We initially used the classic compressibility formula |
| 433 |
> |
Isothermal compressibility values calculated using the number density |
| 434 |
> |
(Eq. \ref{eq:BMN}) expression are in good agreement with experimental |
| 435 |
> |
and previous simulation work throughout the 1 - 1000 atm pressure |
| 436 |
> |
regime. Compressibilities computed using the Hull volume, however, |
| 437 |
> |
deviate dramatically from the experimental values at low applied |
| 438 |
> |
pressures. The reason for this deviation is quite simple; at low |
| 439 |
> |
applied pressures, the liquid is in equilibrium with a vapor phase, |
| 440 |
> |
and it is entirely possible for one (or a few) molecules to drift away |
| 441 |
> |
from the liquid cluster (see Fig. \ref{fig:coneOfShame}). At low |
| 442 |
> |
pressures, the restoring forces on the facets are very gentle, and |
| 443 |
> |
this means that the hulls often take on relatively distorted |
| 444 |
> |
geometries which include large volumes of empty space. |
| 445 |
|
|
| 446 |
< |
\begin{equation} |
| 447 |
< |
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T} |
| 448 |
< |
\end{equation} |
| 446 |
> |
\begin{figure} |
| 447 |
> |
\includegraphics[width=\linewidth]{flytest2} |
| 448 |
> |
\caption{At low pressures, the liquid is in equilibrium with the vapor |
| 449 |
> |
phase, and isolated molecules can detach from the liquid droplet. |
| 450 |
> |
This is expected behavior, but the reported volume of the convex |
| 451 |
> |
hull includes large regions of empty space. For this reason, |
| 452 |
> |
compressibilities should be computed using local number densities |
| 453 |
> |
rather than hull volumes.} |
| 454 |
> |
\label{fig:coneOfShame} |
| 455 |
> |
\end{figure} |
| 456 |
|
|
| 457 |
< |
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. |
| 457 |
> |
At higher pressures, the equilibrium favors the liquid phase, and the |
| 458 |
> |
hull geometries are much more compact. Because of the liquid-vapor |
| 459 |
> |
effect on the convex hull, the regional number density approach |
| 460 |
> |
(Eq. \ref{eq:BMN}) provides more reliable estimates of the bulk |
| 461 |
> |
modulus. |
| 462 |
|
|
| 463 |
+ |
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. |
| 464 |
+ |
|
| 465 |
+ |
Regardless of the difficulty in obtaining accurate hull |
| 466 |
+ |
volumes at low temperature and pressures, the Langevin Hull NPT method |
| 467 |
+ |
provides reasonable isothermal compressibility values for water |
| 468 |
+ |
through a large range of pressures. |
| 469 |
+ |
|
| 470 |
|
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 |
| 471 |
|
|
| 472 |
|
\begin{equation} |
| 475 |
|
|
| 476 |
|
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. |
| 477 |
|
|
| 403 |
– |
In order to calculate the isothermal compressibility without being hindered by hull volume issues, we adapted the classic compressibility formula so that the compressibility could be calculated using information about the local density instead of the volume of the convex hull. We calculated the $g_{OO}(r)$ for a 1 nanosecond simulation of a cluster of 1372 SPC/E water molecules and spherically integrated the function over the bounds 0 to $r'$. In all cases, the value of $r'$ was 17.26216 $\AA$. The value of the total integral between these bounds is essentially the number (N) of molecules within volume $\frac{4}{3}\pi r'^{3}$ at a given pressure. To yield an actual molecule count, N must be scaled by an ideal density. However, even in the absence of an ideal density, we can use the relationship $\rho = \frac{N}{V}$ to rewrite the isothermal compressibility formula as |
| 478 |
|
|
| 405 |
– |
\begin{equation} |
| 406 |
– |
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T} |
| 407 |
– |
\end{equation} |
| 408 |
– |
|
| 409 |
– |
Isothermal compressibility values calculated using this modified expression are in good agreement with the reference values throughout the 1 - 1000 atm pressure regime. Regardless of the difficulty in obtaining accurate hull volumes at low temperature and pressures, the Langevin Hull NPT method provides reasonable isothermal compressibility values for water through a large range of pressures. |
| 410 |
– |
|
| 479 |
|
\subsection{Molecular orientation distribution at cluster boundary} |
| 480 |
|
|
| 481 |
|
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. |
