1 |
gezelter |
3640 |
\documentclass[11pt]{article} |
2 |
|
|
\usepackage{amsmath} |
3 |
|
|
\usepackage{amssymb} |
4 |
|
|
\usepackage{setspace} |
5 |
|
|
\usepackage{endfloat} |
6 |
|
|
\usepackage{caption} |
7 |
|
|
\usepackage{graphicx} |
8 |
|
|
\usepackage{multirow} |
9 |
|
|
\usepackage[square, comma, sort&compress]{natbib} |
10 |
|
|
\usepackage{url} |
11 |
|
|
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm |
12 |
|
|
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight |
13 |
|
|
9.0in \textwidth 6.5in \brokenpenalty=10000 |
14 |
|
|
|
15 |
|
|
% double space list of tables and figures |
16 |
|
|
%\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}} |
17 |
|
|
\setlength{\abovecaptionskip}{20 pt} |
18 |
|
|
\setlength{\belowcaptionskip}{30 pt} |
19 |
|
|
|
20 |
|
|
\bibpunct{[}{]}{,}{s}{}{;} |
21 |
gezelter |
3667 |
\bibliographystyle{achemso} |
22 |
gezelter |
3640 |
|
23 |
|
|
\begin{document} |
24 |
|
|
|
25 |
|
|
\title{The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems} |
26 |
|
|
|
27 |
kstocke1 |
3644 |
\author{Charles F. Vardeman II, Kelsey M. Stocker, and J. Daniel |
28 |
gezelter |
3640 |
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ |
29 |
|
|
Department of Chemistry and Biochemistry,\\ |
30 |
|
|
University of Notre Dame\\ |
31 |
|
|
Notre Dame, Indiana 46556} |
32 |
|
|
|
33 |
|
|
\date{\today} |
34 |
|
|
|
35 |
|
|
\maketitle |
36 |
|
|
|
37 |
|
|
\begin{doublespace} |
38 |
|
|
|
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 |
gezelter |
3665 |
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 |
gezelter |
3652 |
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 |
gezelter |
3665 |
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 |
gezelter |
3640 |
\end{abstract} |
56 |
|
|
|
57 |
|
|
\newpage |
58 |
|
|
|
59 |
|
|
%\narrowtext |
60 |
|
|
|
61 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
62 |
|
|
% BODY OF TEXT |
63 |
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
64 |
|
|
|
65 |
|
|
|
66 |
|
|
\section{Introduction} |
67 |
|
|
|
68 |
gezelter |
3641 |
The most common molecular dynamics methods for sampling configurations |
69 |
gezelter |
3667 |
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{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 |
gezelter |
3665 |
coordinate transformation to adjust the box volume. As long as the |
85 |
gezelter |
3667 |
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 |
gezelter |
3652 |
|
90 |
gezelter |
3665 |
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 |
gezelter |
3652 |
interfacial environments. In these cases, the affine transform of |
96 |
|
|
atomic coordinates will either cause numerical instability when the |
97 |
gezelter |
3665 |
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 |
gezelter |
3652 |
|
101 |
gezelter |
3640 |
\begin{figure} |
102 |
gezelter |
3641 |
\includegraphics[width=\linewidth]{AffineScale2} |
103 |
gezelter |
3667 |
\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 |
gezelter |
3640 |
\label{affineScale} |
111 |
|
|
\end{figure} |
112 |
|
|
|
113 |
gezelter |
3653 |
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 |
gezelter |
3665 |
volume requires either effective solute concentrations that are much |
117 |
gezelter |
3653 |
higher than desirable, or unreasonable system sizes to avoid this |
118 |
gezelter |
3665 |
effect. For example, calculations using typical hydration shells |
119 |
gezelter |
3653 |
solvating a protein under periodic boundary conditions are quite |
120 |
|
|
expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL |
121 |
|
|
SIMULATIONS] |
122 |
gezelter |
3640 |
|
123 |
gezelter |
3665 |
\subsection*{Boundary Methods} |
124 |
gezelter |
3667 |
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 |
gezelter |
3640 |
|
136 |
gezelter |
3653 |
The electrostatic and dispersive behavior near the boundary has long |
137 |
gezelter |
3665 |
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 |
gezelter |
3640 |
|
157 |
gezelter |
3665 |
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 |
gezelter |
3640 |
|
163 |
gezelter |
3665 |
\subsection*{Restraining Potentials} |
164 |
gezelter |
3653 |
Restraining {\it potentials} introduce repulsive potentials at the |
165 |
|
|
surface of a sphere or other geometry. The solute and any explicit |
166 |
gezelter |
3665 |
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 |
gezelter |
3653 |
|
173 |
gezelter |
3665 |
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 |
gezelter |
3653 |
identical to a fixed-volume restraining potential. |
185 |
|
|
|
186 |
gezelter |
3665 |
\subsection*{Hull methods} |
187 |
gezelter |
3653 |
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 |
|
|
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 |
gezelter |
3667 |
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 |
gezelter |
3653 |
|
213 |
gezelter |
3640 |
\section{Methodology} |
214 |
gezelter |
3653 |
\label{sec:meth} |
215 |
gezelter |
3640 |
|
216 |
gezelter |
3665 |
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 |
gezelter |
3640 |
|
229 |
gezelter |
3652 |
The convex hull is the set of facets that have {\it no concave |
230 |
gezelter |
3665 |
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 |
gezelter |
3660 |
with the atoms on the edge and not with atoms interior to the |
238 |
|
|
simulation. |
239 |
gezelter |
3640 |
|
240 |
gezelter |
3662 |
\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 |
gezelter |
3667 |
computed dynamically at each time step, and molecules can move |
245 |
|
|
between the interior (Newtonian) region and the Langevin hull.} |
246 |
gezelter |
3662 |
\label{fig:hullSample} |
247 |
|
|
\end{figure} |
248 |
|
|
|
249 |
gezelter |
3665 |
Atomic sites in the interior of the simulation move under standard |
250 |
gezelter |
3660 |
Newtonian dynamics, |
251 |
gezelter |
3640 |
\begin{equation} |
252 |
gezelter |
3652 |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
253 |
|
|
\label{eq:Newton} |
254 |
gezelter |
3640 |
\end{equation} |
255 |
gezelter |
3652 |
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 |
gezelter |
3640 |
\begin{equation} |
262 |
gezelter |
3652 |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
263 |
gezelter |
3640 |
\end{equation} |
264 |
|
|
|
265 |
gezelter |
3665 |
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 |
gezelter |
3640 |
\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 |
gezelter |
3652 |
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 |
gezelter |
3660 |
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 |
gezelter |
3640 |
\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 |
gezelter |
3665 |
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 |
gezelter |
3652 |
\begin{equation} |
292 |
|
|
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
293 |
|
|
\end{equation} |
294 |
gezelter |
3660 |
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 |
gezelter |
3640 |
\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 |
gezelter |
3652 |
\Xi_f(t)\delta(t-t^\prime). |
303 |
|
|
\label{eq:randomForce} |
304 |
gezelter |
3640 |
\end{eqnarray} |
305 |
|
|
|
306 |
gezelter |
3665 |
Once the resistance tensor is known for a given facet, a stochastic |
307 |
gezelter |
3660 |
vector that has the properties in Eq. (\ref{eq:randomForce}) can be |
308 |
gezelter |
3665 |
calculated efficiently by carrying out a Cholesky decomposition to |
309 |
|
|
obtain the square root matrix of the resistance tensor, |
310 |
gezelter |
3652 |
\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 |
gezelter |
3640 |
|
326 |
gezelter |
3660 |
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 |
gezelter |
3653 |
tensor that includes translational and rotational drag as well as |
329 |
gezelter |
3660 |
translational-rotational coupling. The computation of resistance |
330 |
gezelter |
3653 |
tensors for rigid bodies has been detailed |
331 |
gezelter |
3663 |
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun:2008fk} |
332 |
gezelter |
3653 |
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 |
gezelter |
3665 |
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 |
gezelter |
3653 |
\begin{equation} |
340 |
gezelter |
3665 |
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 |
gezelter |
3653 |
\end{equation} |
343 |
gezelter |
3665 |
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 |
gezelter |
3653 |
|
349 |
|
|
\begin{figure} |
350 |
|
|
\includegraphics[width=\linewidth]{hydro} |
351 |
gezelter |
3660 |
\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 |
gezelter |
3653 |
easily computed using half the cross product of two of the edges.} |
357 |
|
|
\label{hydro} |
358 |
|
|
\end{figure} |
359 |
|
|
|
360 |
gezelter |
3665 |
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 |
gezelter |
3653 |
\begin{equation} |
364 |
|
|
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
365 |
|
|
\end{equation} |
366 |
gezelter |
3667 |
Note that this treatment ignores rotations (and |
367 |
gezelter |
3660 |
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 |
gezelter |
3652 |
We have implemented this method by extending the Langevin dynamics |
372 |
gezelter |
3665 |
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 |
gezelter |
3667 |
\item Delaunay triangulation is carried out using the current atomic |
377 |
gezelter |
3665 |
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 |
gezelter |
3652 |
|
400 |
gezelter |
3640 |
\section{Tests \& Applications} |
401 |
gezelter |
3653 |
\label{sec:tests} |
402 |
gezelter |
3640 |
|
403 |
gezelter |
3663 |
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 |
gezelter |
3665 |
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 |
gezelter |
3660 |
\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 |
gezelter |
3665 |
\frac{N}{V}$, is uniform within some region of the point cloud. The |
421 |
gezelter |
3660 |
compressibility can then be expressed in terms of the average number |
422 |
|
|
of particles in that region, |
423 |
|
|
\begin{equation} |
424 |
gezelter |
3665 |
\kappa_{T} = -\frac{1}{N} \left ( \frac{\partial N}{\partial P} \right |
425 |
gezelter |
3660 |
)_{T} |
426 |
|
|
\label{eq:BMN} |
427 |
|
|
\end{equation} |
428 |
gezelter |
3663 |
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 |
gezelter |
3665 |
cluster can be used to compute the compressibility. |
433 |
gezelter |
3660 |
|
434 |
gezelter |
3665 |
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 |
gezelter |
3660 |
|
439 |
gezelter |
3663 |
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 |
gezelter |
3667 |
\subsection{Compressibility of gold nanoparticles} |
451 |
gezelter |
3640 |
|
452 |
gezelter |
3665 |
The compressibility is well-known for gold, and it provides a good first |
453 |
gezelter |
3663 |
test of how the method compares to other similar methods. |
454 |
|
|
|
455 |
gezelter |
3640 |
\begin{figure} |
456 |
gezelter |
3665 |
\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 |
gezelter |
3640 |
\label{pressureResponse} |
465 |
|
|
\end{figure} |
466 |
|
|
|
467 |
|
|
\begin{equation} |
468 |
|
|
\kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial |
469 |
|
|
P}\right) |
470 |
|
|
\end{equation} |
471 |
|
|
|
472 |
|
|
\begin{figure} |
473 |
|
|
\includegraphics[width=\linewidth]{compress_tb} |
474 |
|
|
\caption{Isothermal Compressibility (18 \AA gold nanoparticle)} |
475 |
|
|
\label{temperatureResponse} |
476 |
|
|
\end{figure} |
477 |
|
|
|
478 |
|
|
\subsection{Compressibility of SPC/E water clusters} |
479 |
|
|
|
480 |
gezelter |
3660 |
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 |
kstocke1 |
3649 |
|
491 |
gezelter |
3640 |
\begin{figure} |
492 |
gezelter |
3659 |
\includegraphics[width=\linewidth]{new_isothermalN} |
493 |
kstocke1 |
3649 |
\caption{Compressibility of SPC/E water} |
494 |
gezelter |
3660 |
\label{fig:compWater} |
495 |
gezelter |
3640 |
\end{figure} |
496 |
|
|
|
497 |
gezelter |
3660 |
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 |
kstocke1 |
3649 |
|
510 |
gezelter |
3660 |
\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 |
gezelter |
3665 |
This is expected behavior, but the volume of the convex hull |
515 |
|
|
includes large regions of empty space. For this reason, |
516 |
gezelter |
3662 |
compressibilities are computed using local number densities rather |
517 |
|
|
than hull volumes.} |
518 |
gezelter |
3660 |
\label{fig:coneOfShame} |
519 |
|
|
\end{figure} |
520 |
kstocke1 |
3649 |
|
521 |
gezelter |
3665 |
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 |
gezelter |
3667 |
compressibility. |
526 |
kstocke1 |
3649 |
|
527 |
gezelter |
3665 |
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 |
kstocke1 |
3649 |
\begin{equation} |
534 |
gezelter |
3665 |
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle |
535 |
|
|
V \right \rangle ^{2}}{V \, k_{B} \, T}, |
536 |
kstocke1 |
3649 |
\end{equation} |
537 |
gezelter |
3665 |
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 |
gezelter |
3667 |
magnitude larger than the reference values. However, compressibility |
547 |
gezelter |
3665 |
calculation that relies on the hull volume will suffer these effects. |
548 |
|
|
WE NEED MORE HERE. |
549 |
kstocke1 |
3649 |
|
550 |
|
|
\subsection{Molecular orientation distribution at cluster boundary} |
551 |
|
|
|
552 |
gezelter |
3665 |
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 |
gezelter |
3667 |
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 |
kstocke1 |
3649 |
|
562 |
gezelter |
3667 |
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 |
kstocke1 |
3649 |
|
574 |
gezelter |
3667 |
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 |
kstocke1 |
3649 |
\begin{equation} |
579 |
gezelter |
3640 |
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|} |
580 |
|
|
\end{equation} |
581 |
gezelter |
3667 |
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 |
gezelter |
3640 |
\begin{figure} |
586 |
kstocke1 |
3649 |
\includegraphics[width=\linewidth]{g_r_theta} |
587 |
gezelter |
3667 |
\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 |
kstocke1 |
3649 |
\end{figure} |
595 |
|
|
|
596 |
gezelter |
3667 |
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 |
kstocke1 |
3649 |
\begin{figure} |
600 |
gezelter |
3640 |
\includegraphics[width=\linewidth]{pAngle} |
601 |
gezelter |
3667 |
\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 |
gezelter |
3640 |
\end{figure} |
609 |
|
|
|
610 |
gezelter |
3667 |
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 |
gezelter |
3640 |
|
617 |
gezelter |
3667 |
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 |
kstocke1 |
3649 |
|
632 |
gezelter |
3667 |
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 |
kstocke1 |
3649 |
|
638 |
gezelter |
3640 |
\subsection{Heterogeneous nanoparticle / water mixtures} |
639 |
|
|
|
640 |
gezelter |
3665 |
\section{Discussion} |
641 |
|
|
\label{sec:discussion} |
642 |
gezelter |
3640 |
|
643 |
gezelter |
3667 |
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 |
gezelter |
3663 |
\section*{Appendix A: Computing Convex Hulls on Parallel Computers} |
661 |
gezelter |
3640 |
|
662 |
gezelter |
3666 |
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 |
gezelter |
3667 |
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 |
gezelter |
3666 |
\begin{enumerate} |
676 |
|
|
\item Each processor computes the convex hull for its own atomic sites |
677 |
|
|
(dashed lines 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 |
gezelter |
3667 |
\item Each processor computes the global convex hull (solid line 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 |
gezelter |
3666 |
\end{enumerate} |
687 |
|
|
|
688 |
|
|
\begin{figure} |
689 |
|
|
\begin{centering} |
690 |
|
|
\includegraphics[width=3in]{parallel} |
691 |
|
|
\caption{When the sites are distributed among many nodes for parallel |
692 |
|
|
computation, the processors first compute the convex hulls for their |
693 |
|
|
own sites (dashed lines in upper panel). The positions of the sites |
694 |
|
|
that make up the convex hulls are then communicated to all |
695 |
|
|
processors. The convex hull of the system (solid line in lower |
696 |
|
|
panel) is the convex hull of the points on the hulls for all |
697 |
|
|
processors. } |
698 |
|
|
\end{centering} |
699 |
|
|
\label{fig:parallel} |
700 |
|
|
\end{figure} |
701 |
|
|
|
702 |
|
|
The individual hull operations scale with |
703 |
gezelter |
3667 |
$\mathcal{O}(\frac{n}{p}\log\frac{n}{p})$ where $n$ is the total |
704 |
|
|
number of sites, and $p$ is the number of processors. These local |
705 |
|
|
hull operations create a set of $p$ hulls each with approximately |
706 |
|
|
$\frac{n}{3pr}$ sites (for a cluster of radius $r$). The worst-case |
707 |
|
|
communication cost for using a ``gather'' operation to distribute this |
708 |
|
|
information to all processors is $\mathcal{O}( \alpha (p-1) + \frac{n |
709 |
|
|
\beta (p-1)}{3 r p^2})$, while the final computation of the system |
710 |
|
|
hull scales as $\mathcal{O}(\frac{n}{3r}\log\frac{n}{3r})$. |
711 |
gezelter |
3666 |
|
712 |
gezelter |
3667 |
For a large number of atoms on a moderately parallel machine, the |
713 |
|
|
total costs are dominated by the computations of the individual hulls, |
714 |
|
|
and communication of these hulls to so the Langevin hull sees roughly |
715 |
|
|
linear speed-up with increasing processor counts. |
716 |
|
|
|
717 |
gezelter |
3663 |
\section*{Acknowledgments} |
718 |
gezelter |
3640 |
Support for this project was provided by the |
719 |
|
|
National Science Foundation under grant CHE-0848243. Computational |
720 |
|
|
time was provided by the Center for Research Computing (CRC) at the |
721 |
|
|
University of Notre Dame. |
722 |
|
|
|
723 |
|
|
\newpage |
724 |
|
|
|
725 |
|
|
\bibliography{langevinHull} |
726 |
|
|
|
727 |
|
|
\end{doublespace} |
728 |
|
|
\end{document} |