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 |
> |
\begin{centering} |
690 |
> |
\includegraphics[width=\linewidth]{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 left panel). The positions of the sites |
694 |
> |
that make up the convex hulls are then communicated to all |
695 |
> |
processors (middle panel). The convex hull of the system (solid line in right panel) is the convex hull of the points on the hulls for all |
696 |
> |
processors.} |
697 |
> |
\label{fig:parallel} |
698 |
> |
\end{centering} |
699 |
> |
\label{fig:parallel} |
700 |
|
\end{figure} |
701 |
|
|
702 |
< |
\begin{equation} |
703 |
< |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1} |
704 |
< |
\end{equation} |
702 |
> |
The individual hull operations scale with |
703 |
> |
$\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 |
|
|
712 |
< |
\begin{equation} |
713 |
< |
T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I + |
714 |
< |
\frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right) |
715 |
< |
\end{equation} |
712 |
> |
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 |
< |
\section{Appendix B: Computing Convex Hulls on Parallel Computers} |
245 |
< |
|
246 |
< |
\section{Acknowledgments} |
717 |
> |
\section*{Acknowledgments} |
718 |
|
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 |