| 18 |
|
\setlength{\belowcaptionskip}{30 pt} |
| 19 |
|
|
| 20 |
|
\bibpunct{[}{]}{,}{s}{}{;} |
| 21 |
< |
\bibliographystyle{aip} |
| 21 |
> |
\bibliographystyle{achemso} |
| 22 |
|
|
| 23 |
|
\begin{document} |
| 24 |
|
|
| 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'', |
| 45 |
< |
performs better than traditional affine transform methods for |
| 46 |
< |
systems containing heterogeneous mixtures of materials with |
| 47 |
< |
different compressibilities. It does not suffer from the edge |
| 48 |
< |
effects of boundary potential methods, and allows realistic |
| 49 |
< |
treatment of both external pressure and thermal conductivity to an |
| 50 |
< |
implicit solvents. We apply this method to several different |
| 51 |
< |
systems including bare nano-particles, nano-particles in explicit |
| 52 |
< |
solvent, as well as clusters of liquid water and ice. The predicted |
| 53 |
< |
mechanical and thermal properties of these systems are in good |
| 54 |
< |
agreement with experimental data. |
| 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 |
> |
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 |
| 69 |
> |
from an isobaric-isothermal (NPT) ensemble maintain a target pressure |
| 70 |
> |
in a simulation by coupling the volume of the system to a {\it |
| 71 |
> |
barostat}, which is an extra degree of freedom propagated along with |
| 72 |
> |
the particle coordinates. These methods require periodic boundary |
| 73 |
> |
conditions, because when the instantaneous pressure in the system |
| 74 |
> |
differs from the target pressure, the volume is reduced or expanded |
| 75 |
> |
using {\it affine transforms} of the system geometry. An affine |
| 76 |
> |
transform scales the size and shape of the periodic box as well as the |
| 77 |
> |
particle positions within the box (but not the sizes of the |
| 78 |
|
particles). The most common constant pressure methods, including the |
| 79 |
< |
Melchionna modification\cite{melchionna93} to the |
| 80 |
< |
Nos\'e-Hoover-Andersen equations of motion, the Berendsen pressure |
| 81 |
< |
bath, and the Langevin Piston, all utilize coordinate transformation |
| 82 |
< |
to adjust the box volume. |
| 79 |
> |
Melchionna modification\cite{Melchionna1993} to the |
| 80 |
> |
Nos\'e-Hoover-Andersen equations of |
| 81 |
> |
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen |
| 82 |
> |
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin |
| 83 |
> |
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize scaled |
| 84 |
> |
coordinate transformation to adjust the box volume. As long as the |
| 85 |
> |
material in the simulation box has a relatively uniform |
| 86 |
> |
compressibility, the standard affine transform approach provides an |
| 87 |
> |
excellent way of adjusting the volume of the system and applying |
| 88 |
> |
pressure directly via the interactions between atomic sites. |
| 89 |
|
|
| 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 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} |
| 103 |
< |
\caption{Affine Scaling constant pressure methods use box-length |
| 104 |
< |
scaling to adjust the volume to adjust to under- or over-pressure |
| 105 |
< |
conditions. In a system with a uniform compressibility (e.g. bulk |
| 106 |
< |
fluids) these methods can work well. In systems containing |
| 107 |
< |
heterogeneous mixtures, the affine scaling moves required to adjust |
| 108 |
< |
the pressure in the high-compressibility regions can cause molecules |
| 109 |
< |
in low compressibility regions to collide.} |
| 103 |
> |
\caption{Affine scaling methods use box-length scaling to adjust the |
| 104 |
> |
volume to adjust to under- or over-pressure conditions. In a system |
| 105 |
> |
with a uniform compressibility (e.g. bulk fluids) these methods can |
| 106 |
> |
work well. In systems containing heterogeneous mixtures, the affine |
| 107 |
> |
scaling moves required to adjust the pressure in the |
| 108 |
> |
high-compressibility regions can cause molecules in low |
| 109 |
> |
compressibility regions to collide.} |
| 110 |
|
\label{affineScale} |
| 111 |
|
\end{figure} |
| 112 |
|
|
| 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 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 |
| 119 |
+ |
solvating a protein under periodic boundary conditions are quite |
| 120 |
+ |
expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL |
| 121 |
+ |
SIMULATIONS] |
| 122 |
|
|
| 123 |
< |
Heterogeneous mixtures of materials with different compressibilities? |
| 123 |
> |
\subsection*{Boundary Methods} |
| 124 |
> |
There have been a number of approaches to handle simulations of |
| 125 |
> |
explicitly non-periodic systems that focus on constant or |
| 126 |
> |
nearly-constant {\it volume} conditions while maintaining bulk-like |
| 127 |
> |
behavior. Berkowitz and McCammon introduced a stochastic (Langevin) |
| 128 |
> |
boundary layer inside a region of fixed molecules which effectively |
| 129 |
> |
enforces constant temperature and volume (NVT) |
| 130 |
> |
conditions.\cite{Berkowitz1982} In this approach, the stochastic and |
| 131 |
> |
fixed regions were defined relative to a central atom. Brooks and |
| 132 |
> |
Karplus extended this method to include deformable stochastic |
| 133 |
> |
boundaries.\cite{iii:6312} The stochastic boundary approach has been |
| 134 |
> |
used widely for protein simulations. [CITATIONS NEEDED] |
| 135 |
|
|
| 136 |
< |
Explicitly non-periodic systems |
| 136 |
> |
The electrostatic and dispersive behavior near the boundary has long |
| 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 |
< |
Elastic Bag |
| 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 |
< |
Spherical Boundary approaches |
| 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 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 |
< |
\section{Methodology} |
| 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 |
< |
A new method which uses a constant pressure and temperature bath that |
| 187 |
< |
interacts with the objects that are currently at the edge of the |
| 188 |
< |
system. |
| 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} |
| 190 |
> |
This method is based on standard Langevin dynamics, but the Brownian |
| 191 |
> |
or random forces are allowed to act only on peripheral atoms and exert |
| 192 |
> |
force in a direction that is inward-facing relative to the facets of a |
| 193 |
> |
closed bounding surface. The statistical distribution of the random |
| 194 |
> |
forces are uniquely tied to the pressure in the external reservoir, so |
| 195 |
> |
the method can be shown to sample the isobaric-isothermal ensemble. |
| 196 |
> |
Kohanoff {\it et al.} used a Delaunay tessellation to generate a |
| 197 |
> |
bounding surface surrounding the outermost atoms in the simulated |
| 198 |
> |
system. This is not the only possible triangulated outer surface, but |
| 199 |
> |
guarantees that all of the random forces point inward towards the |
| 200 |
> |
cluster. |
| 201 |
|
|
| 202 |
< |
Novel features: No a priori geometry is defined, No affine transforms, |
| 203 |
< |
No fictitious particles, No bounding potentials. |
| 204 |
< |
|
| 205 |
< |
Simulation starts as a collection of atomic locations in 3D (a point |
| 206 |
< |
cloud). |
| 202 |
> |
In the following sections, we extend and generalize the approach of |
| 203 |
> |
Kohanoff, Caro, and Finnis. The new method, which we are calling the |
| 204 |
> |
``Langevin Hull'' applies the external pressure, Langevin drag, and |
| 205 |
> |
random forces on the {\it facets of the hull} instead of the atomic |
| 206 |
> |
sites comprising the vertices of the hull. This allows us to decouple |
| 207 |
> |
the external pressure contribution from the drag and random force. |
| 208 |
> |
The methodology is introduced in section \ref{sec:meth}, tests on |
| 209 |
> |
crystalline nanoparticles, liquid clusters, and heterogeneous mixtures |
| 210 |
> |
are detailed in section \ref{sec:tests}. Section \ref{sec:discussion} |
| 211 |
> |
summarizes our findings. |
| 212 |
|
|
| 213 |
< |
Delaunay triangulation finds all facets between coplanar neighbors. |
| 213 |
> |
\section{Methodology} |
| 214 |
> |
\label{sec:meth} |
| 215 |
|
|
| 216 |
< |
The Convex Hull is the set of facets that have no concave corners at a |
| 217 |
< |
vertex. |
| 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 |
< |
Molecules on the convex hull are dynamic. As they re-enter the |
| 230 |
< |
cluster, all interactions with the external bath are removed.The |
| 231 |
< |
external bath applies pressure to the facets of the convex hull in |
| 232 |
< |
direct proportion to the area of the facet.Thermal coupling depends on |
| 233 |
< |
the solvent temperature, friction and the size and shape of each |
| 234 |
< |
facet. |
| 229 |
> |
The convex hull is the set of facets that have {\it no concave |
| 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 |
|
|
| 240 |
+ |
\begin{figure} |
| 241 |
+ |
\includegraphics[width=\linewidth]{hullSample} |
| 242 |
+ |
\caption{The external temperature and pressure bath interacts only |
| 243 |
+ |
with those atoms on the convex hull (grey surface). The hull is |
| 244 |
+ |
computed dynamically at each time step, and molecules can move |
| 245 |
+ |
between the interior (Newtonian) region and the Langevin hull.} |
| 246 |
+ |
\label{fig:hullSample} |
| 247 |
+ |
\end{figure} |
| 248 |
+ |
|
| 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 |
| 252 |
> |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
| 253 |
> |
\label{eq:Newton} |
| 254 |
|
\end{equation} |
| 255 |
< |
|
| 255 |
> |
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the |
| 256 |
> |
instantaneous velocity of site $i$ at time $t$, and $U$ is the total |
| 257 |
> |
potential energy. For atoms on the exterior of the cluster |
| 258 |
> |
(i.e. those that occupy one of the vertices of the convex hull), the |
| 259 |
> |
equation of motion is modified with an external force, ${\mathbf |
| 260 |
> |
F}_i^{\mathrm ext}$, |
| 261 |
|
\begin{equation} |
| 262 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext} |
| 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 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 |
| 274 |
|
F}_f^{\mathrm ext} |
| 275 |
|
\end{equation} |
| 276 |
|
|
| 277 |
+ |
The external pressure bath applies a force to the facets of the convex |
| 278 |
+ |
hull in direct proportion to the area of the facet, while the thermal |
| 279 |
+ |
coupling depends on the solvent temperature, viscosity and the size |
| 280 |
+ |
and shape of each facet. The thermal interactions are expressed as a |
| 281 |
+ |
standard Langevin description of the forces, |
| 282 |
|
\begin{equation} |
| 283 |
|
\begin{array}{rclclcl} |
| 284 |
|
{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\ |
| 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 |
< |
|
| 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} |
| 294 |
> |
and $\Xi_f(t)$ is an approximate ($3 \times 3$) resistance tensor that |
| 295 |
> |
depends on the geometry and surface area of facet $f$ and the |
| 296 |
> |
viscosity of the fluid. The resistance tensor is related to the |
| 297 |
> |
fluctuations of the random force, $\mathbf{R}(t)$, by the |
| 298 |
> |
fluctuation-dissipation theorem, |
| 299 |
|
\begin{eqnarray} |
| 150 |
– |
A_f & = & \text{area of facet}\ f \\ |
| 151 |
– |
\hat{n}_f & = & \text{facet normal} \\ |
| 152 |
– |
P & = & \text{external pressure} |
| 153 |
– |
\end{eqnarray} |
| 154 |
– |
|
| 155 |
– |
\begin{eqnarray} |
| 156 |
– |
{\mathbf v}_f(t) & = & \text{velocity of facet} \\ |
| 157 |
– |
& = & \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i \\ |
| 158 |
– |
\Xi_f(t) & = & \text{is a hydrodynamic tensor that depends} \\ |
| 159 |
– |
& & \text{on the geometry and surface area of} \\ |
| 160 |
– |
& & \text{facet} \ f\ \text{and the viscosity of the fluid.} |
| 161 |
– |
\end{eqnarray} |
| 162 |
– |
|
| 163 |
– |
\begin{eqnarray} |
| 300 |
|
\left< {\mathbf R}_f(t) \right> & = & 0 \\ |
| 301 |
|
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\ |
| 302 |
< |
\Xi_f(t)\delta(t-t^\prime) |
| 302 |
> |
\Xi_f(t)\delta(t-t^\prime). |
| 303 |
> |
\label{eq:randomForce} |
| 304 |
|
\end{eqnarray} |
| 305 |
|
|
| 306 |
< |
Implemented in OpenMD.\cite{Meineke:2005gd,openmd} |
| 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 |
> |
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} |
| 313 |
> |
\end{equation} |
| 314 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
| 315 |
> |
vector with the statistics required for the random force can then be |
| 316 |
> |
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which |
| 317 |
> |
has elements chosen from a Gaussian distribution, such that: |
| 318 |
> |
\begin{equation} |
| 319 |
> |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
| 320 |
> |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
| 321 |
> |
\end{equation} |
| 322 |
> |
where $\delta t$ is the timestep in use during the simulation. The |
| 323 |
> |
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
| 324 |
> |
have the correct properties required by Eq. (\ref{eq:randomForce}). |
| 325 |
|
|
| 326 |
< |
\section{Tests \& Applications} |
| 326 |
> |
Our treatment of the resistance tensor is approximate. $\Xi$ for a |
| 327 |
> |
rigid triangular plate would normally be treated as a $6 \times 6$ |
| 328 |
> |
tensor that includes translational and rotational drag as well as |
| 329 |
> |
translational-rotational coupling. The computation of resistance |
| 330 |
> |
tensors for rigid bodies has been detailed |
| 331 |
> |
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun:2008fk} |
| 332 |
> |
but the standard approach involving bead approximations would be |
| 333 |
> |
prohibitively expensive if it were recomputed at each step in a |
| 334 |
> |
molecular dynamics simulation. |
| 335 |
|
|
| 336 |
< |
\subsection{Bulk modulus of gold nanoparticles} |
| 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_{\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_\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]{pressure_tb} |
| 351 |
< |
\caption{Pressure response is rapid (18 \AA gold nanoparticle), target |
| 352 |
< |
pressure = 4 GPa} |
| 353 |
< |
\label{pressureResponse} |
| 350 |
> |
\includegraphics[width=\linewidth]{hydro} |
| 351 |
> |
\caption{The resistance tensor $\Xi$ for a facet comprising sites $i$, |
| 352 |
> |
$j$, and $k$ is constructed using Oseen tensor contributions between |
| 353 |
> |
the centoid of the facet $f$ and each of the sub-facets ($i,f,j$), |
| 354 |
> |
($j,f,k$), and ($k,f,i$). The centroids of the sub-facets are |
| 355 |
> |
located at $1$, $2$, and $3$, and the area of each sub-facet is |
| 356 |
> |
easily computed using half the cross product of two of the edges.} |
| 357 |
> |
\label{hydro} |
| 358 |
|
\end{figure} |
| 359 |
|
|
| 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} |
| 366 |
+ |
Note that this treatment ignores rotations (and |
| 367 |
+ |
translational-rotational coupling) of the facet. In compact systems, |
| 368 |
+ |
the facets stay relatively fixed in orientation between |
| 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} 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 carried out 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{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}. |
| 413 |
+ |
\label{eq:BM} |
| 414 |
+ |
\end{equation} |
| 415 |
+ |
|
| 416 |
+ |
One problem with eliminating periodic boundary conditions and |
| 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 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 |
| 425 |
+ |
)_{T} |
| 426 |
+ |
\label{eq:BMN} |
| 427 |
+ |
\end{equation} |
| 428 |
+ |
The region we used is a spherical volume of 10 \AA\ radius centered in |
| 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 compressibility. |
| 433 |
+ |
|
| 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 |
| 441 |
+ |
simulations involving point charges, we utilized damped shifted-force |
| 442 |
+ |
(DSF) electrostatics\cite{Fennell06} which is a variant of the Wolf |
| 443 |
+ |
summation\cite{wolf:8254} that has been shown to provide good forces |
| 444 |
+ |
and torques on molecular models for water in a computationally |
| 445 |
+ |
efficient manner.\cite{Fennell06} The damping parameter ($\alpha$) was |
| 446 |
+ |
set to 0.18 \AA$^{-1}$, and the cutoff radius was set to 12 \AA. The |
| 447 |
+ |
Spohr potential was adopted in depicting the interaction between metal |
| 448 |
+ |
atoms and the SPC/E water molecules.\cite{ISI:000167766600035} |
| 449 |
+ |
|
| 450 |
+ |
\subsection{Compressibility of gold nanoparticles} |
| 451 |
+ |
|
| 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 |
+ |
|
| 455 |
|
\begin{figure} |
| 456 |
< |
\includegraphics[width=\linewidth]{temperature_tb} |
| 457 |
< |
\caption{Temperature equilibration depends on surface area and bath |
| 458 |
< |
viscosity. Target Temperature = 300K} |
| 459 |
< |
\label{temperatureResponse} |
| 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} |
| 465 |
|
\end{figure} |
| 466 |
|
|
| 467 |
|
\begin{equation} |
| 477 |
|
|
| 478 |
|
\subsection{Compressibility of SPC/E water clusters} |
| 479 |
|
|
| 480 |
+ |
Prior molecular dynamics simulations on SPC/E water (both in |
| 481 |
+ |
NVT~\cite{Glattli2002} and NPT~\cite{Motakabbir1990, Pi2009} |
| 482 |
+ |
ensembles) have yielded values for the isothermal compressibility that |
| 483 |
+ |
agree well with experiment.\cite{Fine1973} The results of two |
| 484 |
+ |
different approaches for computing the isothermal compressibility from |
| 485 |
+ |
Langevin Hull simulations for pressures between 1 and 6500 atm are |
| 486 |
+ |
shown in Fig. \ref{fig:compWater} along with compressibility values |
| 487 |
+ |
obtained from both other SPC/E simulations and experiment. |
| 488 |
+ |
Compressibility values from all references are for applied pressures |
| 489 |
+ |
within the range 1 - 1000 atm. |
| 490 |
+ |
|
| 491 |
|
\begin{figure} |
| 492 |
< |
\includegraphics[width=\linewidth]{g_r_theta} |
| 493 |
< |
\caption{Definition of coordinates} |
| 494 |
< |
\label{coords} |
| 492 |
> |
\includegraphics[width=\linewidth]{new_isothermalN} |
| 493 |
> |
\caption{Compressibility of SPC/E water} |
| 494 |
> |
\label{fig:compWater} |
| 495 |
|
\end{figure} |
| 496 |
|
|
| 497 |
+ |
Isothermal compressibility values calculated using the number density |
| 498 |
+ |
(Eq. \ref{eq:BMN}) expression are in good agreement with experimental |
| 499 |
+ |
and previous simulation work throughout the 1 - 1000 atm pressure |
| 500 |
+ |
regime. Compressibilities computed using the Hull volume, however, |
| 501 |
+ |
deviate dramatically from the experimental values at low applied |
| 502 |
+ |
pressures. The reason for this deviation is quite simple; at low |
| 503 |
+ |
applied pressures, the liquid is in equilibrium with a vapor phase, |
| 504 |
+ |
and it is entirely possible for one (or a few) molecules to drift away |
| 505 |
+ |
from the liquid cluster (see Fig. \ref{fig:coneOfShame}). At low |
| 506 |
+ |
pressures, the restoring forces on the facets are very gentle, and |
| 507 |
+ |
this means that the hulls often take on relatively distorted |
| 508 |
+ |
geometries which include large volumes of empty space. |
| 509 |
+ |
|
| 510 |
+ |
\begin{figure} |
| 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 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 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 |
+ |
compressibility. |
| 526 |
+ |
|
| 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 |
| 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. However, compressibility |
| 547 |
+ |
calculation that relies on the hull volume will suffer these effects. |
| 548 |
+ |
WE NEED MORE HERE. |
| 549 |
+ |
|
| 550 |
+ |
\subsection{Molecular orientation distribution at cluster boundary} |
| 551 |
+ |
|
| 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. Early non-periodic |
| 555 |
+ |
simulation methods (e.g. hydrophobic boundary potentials) induced |
| 556 |
+ |
spurious orientational correlations deep within the simulated |
| 557 |
+ |
system.\cite{Lee1984,Belch1985} This behavior spawned many methods for |
| 558 |
+ |
fixing and characterizing the effects of artifical boundaries |
| 559 |
+ |
including methods which fix the orientations of a set of edge |
| 560 |
+ |
molecules.\cite{Warshel1978,King1989} |
| 561 |
+ |
|
| 562 |
+ |
As described above, the Langevin Hull does not require that the |
| 563 |
+ |
orientation of molecules be fixed, nor does it utilize an explicitly |
| 564 |
+ |
hydrophobic boundary, orientational constraint or radial constraint. |
| 565 |
+ |
Therefore, the orientational correlations of the molecules in a water |
| 566 |
+ |
cluster are of particular interest in testing this method. Ideally, |
| 567 |
+ |
the water molecules on the surface of the cluster will have enough |
| 568 |
+ |
mobility into and out of the center of the cluster to maintain a |
| 569 |
+ |
bulk-like orientational distribution in the absence of orientational |
| 570 |
+ |
and radial constraints. However, since the number of hydrogen bonding |
| 571 |
+ |
partners available to molecules on the exterior are limited, it is |
| 572 |
+ |
likely that there will be some effective hydrophobicity of the hull. |
| 573 |
+ |
|
| 574 |
+ |
To determine the extent of these effects demonstrated by the Langevin |
| 575 |
+ |
Hull, we examined the orientationations exhibited by SPC/E water in a |
| 576 |
+ |
cluster of 1372 molecules at 300 K and at pressures ranging from 1 - |
| 577 |
+ |
1000 atm. The orientational angle of a water molecule is described |
| 578 |
|
\begin{equation} |
| 579 |
|
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|} |
| 580 |
|
\end{equation} |
| 581 |
< |
|
| 581 |
> |
where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of |
| 582 |
> |
mass and the cluster center of mass and $\vec{\mu}_{i}$ is the vector |
| 583 |
> |
bisecting the H-O-H angle of molecule {\it i} (See |
| 584 |
> |
Fig. \ref{fig:coords}). |
| 585 |
|
\begin{figure} |
| 586 |
< |
\includegraphics[width=\linewidth]{pAngle} |
| 587 |
< |
\caption{SPC/E water clusters: only minor dewetting at the boundary} |
| 588 |
< |
\label{pAngle} |
| 586 |
> |
\includegraphics[width=\linewidth]{g_r_theta} |
| 587 |
> |
\caption{Orientation angle of the water molecules relative to the |
| 588 |
> |
center of the cluster. Bulk-like distributions will result in |
| 589 |
> |
$\langle \cos \theta \rangle$ values close to zero. If the hull |
| 590 |
> |
exhibits an overabundance of externally-oriented oxygen sites the |
| 591 |
> |
average orientation will be negative, while dangling hydrogen sites |
| 592 |
> |
will result in positive average orientations.} |
| 593 |
> |
\label{fig:coords} |
| 594 |
|
\end{figure} |
| 595 |
|
|
| 596 |
+ |
Fig. \ref{fig:pAngle} shows the distribution of $\cos{\theta}$ values |
| 597 |
+ |
for molecules in the interior of the cluster (squares) and for |
| 598 |
+ |
molecules included in the convex hull (circles). |
| 599 |
|
\begin{figure} |
| 600 |
< |
\includegraphics[width=\linewidth]{isothermal} |
| 601 |
< |
\caption{Compressibility of SPC/E water} |
| 602 |
< |
\label{compWater} |
| 600 |
> |
\includegraphics[width=\linewidth]{pAngle} |
| 601 |
> |
\caption{Distribution of $\cos{\theta}$ values for molecules on the |
| 602 |
> |
interior of the cluster (squares) and for those participating in the |
| 603 |
> |
convex hull (circles) at a variety of pressures. The Langevin hull |
| 604 |
> |
exhibits minor dewetting behavior with exposed oxygen sites on the |
| 605 |
> |
hull water molecules. The orientational preference for exposed |
| 606 |
> |
oxygen appears to be independent of applied pressure. } |
| 607 |
> |
\label{fig:pAngle} |
| 608 |
|
\end{figure} |
| 609 |
|
|
| 610 |
+ |
As expected, interior molecules (those not included in the convex |
| 611 |
+ |
hull) maintain a bulk-like structure with a uniform distribution of |
| 612 |
+ |
orientations. Molecules included in the convex hull show a slight |
| 613 |
+ |
preference for values of $\cos{\theta} < 0.$ These values correspond |
| 614 |
+ |
to molecules with oxygen directed toward the exterior of the cluster, |
| 615 |
+ |
forming a dangling hydrogen bond acceptor site. |
| 616 |
+ |
|
| 617 |
+ |
In the absence of an electrostatic contribution from the exterior |
| 618 |
+ |
bath, the orientational distribution of water molecules included in |
| 619 |
+ |
the Langevin Hull will slightly resemble the distribution at a neat |
| 620 |
+ |
water liquid/vapor interface. Previous molecular dynamics simulations |
| 621 |
+ |
of SPC/E water \cite{Taylor1996} have shown that molecules at the |
| 622 |
+ |
liquid/vapor interface favor an orientation where one hydrogen |
| 623 |
+ |
protrudes from the liquid phase. This behavior is demonstrated by |
| 624 |
+ |
experiments \cite{Du1994} \cite{Scatena2001} showing that |
| 625 |
+ |
approximately one-quarter of water molecules at the liquid/vapor |
| 626 |
+ |
interface form dangling hydrogen bonds. The negligible preference |
| 627 |
+ |
shown in these cluster simulations could be removed through the |
| 628 |
+ |
introduction of an implicit solvent model, which would provide the |
| 629 |
+ |
missing electrostatic interactions between the cluster molecules and |
| 630 |
+ |
the surrounding temperature/pressure bath. |
| 631 |
+ |
|
| 632 |
+ |
The orientational preference exhibited by hull molecules in the |
| 633 |
+ |
Langevin hull is significantly weaker than the preference caused by an |
| 634 |
+ |
explicit hydrophobic bounding potential. Additionally, the Langevin |
| 635 |
+ |
Hull does not require that the orientation of any molecules be fixed |
| 636 |
+ |
in order to maintain bulk-like structure, even at the cluster surface. |
| 637 |
+ |
|
| 638 |
|
\subsection{Heterogeneous nanoparticle / water mixtures} |
| 639 |
|
|
| 640 |
+ |
\section{Discussion} |
| 641 |
+ |
\label{sec:discussion} |
| 642 |
|
|
| 643 |
< |
\section{Appendix A: Hydrodynamic tensor for triangular facets} |
| 643 |
> |
The Langevin Hull samples the isobaric-isothermal ensemble for |
| 644 |
> |
non-periodic systems by coupling the system to an bath characterized |
| 645 |
> |
by pressure, temperature, and solvent viscosity. This enables the |
| 646 |
> |
study of heterogeneous systems composed of materials of significantly |
| 647 |
> |
different compressibilities. Because the boundary is dynamically |
| 648 |
> |
determined during the simulation and the molecules interacting with |
| 649 |
> |
the boundary can change, the method and has minimal perturbations on |
| 650 |
> |
the behavior of molecules at the edges of the simulation. Further |
| 651 |
> |
work on this method will involve implicit electrostatics at the |
| 652 |
> |
boundary (which is missing in the current implementation) as well as |
| 653 |
> |
more sophisticated treatments of the surface geometry (alpha |
| 654 |
> |
shapes\cite{EDELSBRUNNER:1994oq,EDELSBRUNNER:1995cj} and Tight |
| 655 |
> |
Cocone\cite{Dey:2003ts}). The non-convex hull geometries are |
| 656 |
> |
significantly more expensive ($\mathcal{O}(N^2)$) than the convex hull |
| 657 |
> |
($\mathcal{O}(N \log N)$), but would enable the use of hull volumes |
| 658 |
> |
directly in computing the compressibility of the sample. |
| 659 |
|
|
| 660 |
+ |
\section*{Appendix A: Computing Convex Hulls on Parallel Computers} |
| 661 |
+ |
|
| 662 |
+ |
In order to use the Langevin Hull for simulations on parallel |
| 663 |
+ |
computers, one of the more difficult tasks is to compute the bounding |
| 664 |
+ |
surface, facets, and resistance tensors when the processors have |
| 665 |
+ |
incomplete information about the entire system's topology. Most |
| 666 |
+ |
parallel decomposition methods assign primary responsibility for the |
| 667 |
+ |
motion of an atomic site to a single processor, and we can exploit |
| 668 |
+ |
this to efficiently compute the convex hull for the entire system. |
| 669 |
+ |
|
| 670 |
+ |
The basic idea involves splitting the point cloud into |
| 671 |
+ |
spatially-overlapping subsets and computing the convex hulls for each |
| 672 |
+ |
of the subsets. The points on the convex hull of the entire system |
| 673 |
+ |
are all present on at least one of the subset hulls. The algorithm |
| 674 |
+ |
works as follows: |
| 675 |
+ |
\begin{enumerate} |
| 676 |
+ |
\item Each processor computes the convex hull for its own atomic sites |
| 677 |
+ |
(left panel in Fig. \ref{fig:parallel}). |
| 678 |
+ |
\item The Hull vertices from each processor are passed out to all of |
| 679 |
+ |
the processors, and each processor assembles a complete list of hull |
| 680 |
+ |
sites (this is much smaller than the original number of points in |
| 681 |
+ |
the point cloud). |
| 682 |
+ |
\item Each processor computes the global convex hull (right panel in |
| 683 |
+ |
Fig. \ref{fig:parallel}) using only those points that are the union |
| 684 |
+ |
of sites gathered from all of the subset hulls. Delaunay |
| 685 |
+ |
triangulation is then done to obtain the facets of the global hull. |
| 686 |
+ |
\end{enumerate} |
| 687 |
+ |
|
| 688 |
|
\begin{figure} |
| 689 |
< |
\includegraphics[width=\linewidth]{hydro} |
| 690 |
< |
\caption{Hydro} |
| 691 |
< |
\label{hydro} |
| 689 |
> |
\includegraphics[width=\linewidth]{parallel} |
| 690 |
> |
\caption{When the sites are distributed among many nodes for parallel |
| 691 |
> |
computation, the processors first compute the convex hulls for their |
| 692 |
> |
own sites (dashed lines in left panel). The positions of the sites |
| 693 |
> |
that make up the subset hulls are then communicated to all |
| 694 |
> |
processors (middle panel). The convex hull of the system (solid line in right panel) is the convex hull of the points on the union of the subset hulls.} |
| 695 |
> |
\label{fig:parallel} |
| 696 |
|
\end{figure} |
| 697 |
|
|
| 698 |
< |
\begin{equation} |
| 699 |
< |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1} |
| 700 |
< |
\end{equation} |
| 698 |
> |
The individual hull operations scale with |
| 699 |
> |
$\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total |
| 700 |
> |
number of sites, and $p$ is the number of processors. These local |
| 701 |
> |
hull operations create a set of $p$ hulls each with approximately |
| 702 |
> |
$\frac{n}{3pr}$ sites (for a cluster of radius $r$). The worst-case |
| 703 |
> |
communication cost for using a ``gather'' operation to distribute this |
| 704 |
> |
information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n |
| 705 |
> |
\beta (p-1)}{3 r p^2})$, while the final computation of the system |
| 706 |
> |
hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$. |
| 707 |
|
|
| 708 |
< |
\begin{equation} |
| 709 |
< |
T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I + |
| 710 |
< |
\frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right) |
| 711 |
< |
\end{equation} |
| 708 |
> |
For a large number of atoms on a moderately parallel machine, the |
| 709 |
> |
total costs are dominated by the computations of the individual hulls, |
| 710 |
> |
and communication of these hulls to so the Langevin hull sees roughly |
| 711 |
> |
linear speed-up with increasing processor counts. |
| 712 |
|
|
| 713 |
< |
\section{Appendix B: Computing Convex Hulls on Parallel Computers} |
| 245 |
< |
|
| 246 |
< |
\section{Acknowledgments} |
| 713 |
> |
\section*{Acknowledgments} |
| 714 |
|
Support for this project was provided by the |
| 715 |
|
National Science Foundation under grant CHE-0848243. Computational |
| 716 |
|
time was provided by the Center for Research Computing (CRC) at the |
