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. |
44 |
> |
external heat bath. This new method, the ``Langevin Hull'', performs |
45 |
> |
better than traditional affine transform methods for systems |
46 |
> |
containing heterogeneous mixtures of materials with different |
47 |
> |
compressibilities. It does not suffer from the edge effects of |
48 |
> |
boundary potential methods, and allows realistic treatment of both |
49 |
> |
external pressure and thermal conductivity to an implicit solvent. |
50 |
> |
We apply this method to several different systems including bare |
51 |
> |
nanoparticles, nanoparticles in an explicit solvent, as well as |
52 |
> |
clusters of liquid water and ice. The predicted mechanical and |
53 |
> |
thermal properties of these systems are in good agreement with |
54 |
> |
experimental data. |
55 |
|
\end{abstract} |
56 |
|
|
57 |
|
\newpage |
76 |
|
well as the scaled particle positions (but not the sizes of the |
77 |
|
particles). The most common constant pressure methods, including the |
78 |
|
Melchionna modification\cite{Melchionna1993} to the |
79 |
< |
Nos\'e-Hoover-Andersen equations of motion, the Berendsen pressure |
80 |
< |
bath, and the Langevin Piston, all utilize coordinate transformation |
81 |
< |
to adjust the box volume. |
79 |
> |
Nos\'e-Hoover-Andersen equations of |
80 |
> |
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen |
81 |
> |
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin |
82 |
> |
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize coordinate |
83 |
> |
transformation to adjust the box volume. As long as the material in |
84 |
> |
the simulation box is essentially a bulk-like liquid which has a |
85 |
> |
relatively uniform compressibility, the standard affine transform |
86 |
> |
approach provides an excellent way of adjusting the volume of the |
87 |
> |
system and applying pressure directly via the interactions between |
88 |
> |
atomic sites. |
89 |
|
|
90 |
+ |
The problem with this approach becomes apparent when the material |
91 |
+ |
being simulated is an inhomogeneous mixture in which portions of the |
92 |
+ |
simulation box are incompressible relative to other portions. |
93 |
+ |
Examples include simulations of metallic nanoparticles in liquid |
94 |
+ |
environments, proteins at interfaces, as well as other multi-phase or |
95 |
+ |
interfacial environments. In these cases, the affine transform of |
96 |
+ |
atomic coordinates will either cause numerical instability when the |
97 |
+ |
sites in the incompressible medium collide with each other, or lead to |
98 |
+ |
inefficient sampling of system volumes if the barostat is set slow |
99 |
+ |
enough to avoid the instabilities in the incompressible region. |
100 |
+ |
|
101 |
|
\begin{figure} |
102 |
|
\includegraphics[width=\linewidth]{AffineScale2} |
103 |
|
\caption{Affine Scaling constant pressure methods use box-length |
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 either requires 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 |
> |
There have been a number of other approaches to explicit |
124 |
> |
non-periodicity that focus on constant or nearly-constant {\it volume} |
125 |
> |
conditions while maintaining bulk-like behavior. Berkowitz and |
126 |
> |
McCammon introduced a stochastic (Langevin) boundary layer inside a |
127 |
> |
region of fixed molecules which effectively enforces constant |
128 |
> |
temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this |
129 |
> |
approach, the stochastic and fixed regions were defined relative to a |
130 |
> |
central atom. Brooks and Karplus extended this method to include |
131 |
> |
deformable stochastic boundaries.\cite{iii:6312} The stochastic |
132 |
> |
boundary approach has been used widely for protein |
133 |
> |
simulations. [CITATIONS NEEDED] |
134 |
|
|
135 |
< |
Explicitly non-periodic systems |
135 |
> |
The electrostatic and dispersive behavior near the boundary has long |
136 |
> |
been a cause for concern. King and Warshel introduced a surface |
137 |
> |
constrained all-atom solvent (SCAAS) which included polarization |
138 |
> |
effects of a fixed spherical boundary to mimic bulk-like behavior |
139 |
> |
without periodic boundaries.\cite{king:3647} In the SCAAS model, a |
140 |
> |
layer of fixed solvent molecules surrounds the solute and any explicit |
141 |
> |
solvent, and this in turn is surrounded by a continuum dielectric. |
142 |
> |
MORE HERE. WHAT DID THEY FIND? |
143 |
|
|
144 |
< |
Elastic Bag |
144 |
> |
Beglov and Roux developed a boundary model in which the hard sphere |
145 |
> |
boundary has a radius that varies with the instantaneous configuration |
146 |
> |
of the solute (and solvent) molecules.\cite{beglov:9050} This model |
147 |
> |
contains a clear pressure and surface tension contribution to the free |
148 |
> |
energy which XXX. |
149 |
|
|
150 |
< |
Spherical Boundary approaches |
150 |
> |
Restraining {\it potentials} introduce repulsive potentials at the |
151 |
> |
surface of a sphere or other geometry. The solute and any explicit |
152 |
> |
solvent are therefore restrained inside this potential. Often the |
153 |
> |
potentials include a weak short-range attraction to maintain the |
154 |
> |
correct density at the boundary. Beglov and Roux have also introduced |
155 |
> |
a restraining boundary potential which relaxes dynamically depending |
156 |
> |
on the solute geometry and the force the explicit system exerts on the |
157 |
> |
shell.\cite{Beglov:1995fk} |
158 |
|
|
159 |
< |
\section{Methodology} |
159 |
> |
Recently, Krilov {\it et al.} introduced a flexible boundary model |
160 |
> |
that uses a Lennard-Jones potential between the solvent molecules and |
161 |
> |
a boundary which is determined dynamically from the position of the |
162 |
> |
nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This approach allows |
163 |
> |
the confining potential to prevent solvent molecules from migrating |
164 |
> |
too far from the solute surface, while providing a weak attractive |
165 |
> |
force pulling the solvent molecules towards a fictitious bulk solvent. |
166 |
> |
Although this approach is appealing and has physical motivation, |
167 |
> |
nanoparticles do not deform far from their original geometries even at |
168 |
> |
temperatures which vaporize the nearby solvent. For the systems like |
169 |
> |
the one described, the flexible boundary model will be nearly |
170 |
> |
identical to a fixed-volume restraining potential. |
171 |
|
|
172 |
< |
A new method which uses a constant pressure and temperature bath that |
173 |
< |
interacts with the objects that are currently at the edge of the |
174 |
< |
system. |
172 |
> |
The approach of Kohanoff, Caro, and Finnis is the most promising of |
173 |
> |
the methods for introducing both constant pressure and temperature |
174 |
> |
into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru} |
175 |
> |
This method is based on standard Langevin dynamics, but the Brownian |
176 |
> |
or random forces are allowed to act only on peripheral atoms and exert |
177 |
> |
force in a direction that is inward-facing relative to the facets of a |
178 |
> |
closed bounding surface. The statistical distribution of the random |
179 |
> |
forces are uniquely tied to the pressure in the external reservoir, so |
180 |
> |
the method can be shown to sample the isobaric-isothermal ensemble. |
181 |
> |
Kohanoff {\it et al.} used a Delaunay tessellation to generate a |
182 |
> |
bounding surface surrounding the outermost atoms in the simulated |
183 |
> |
system. This is not the only possible triangulated outer surface, but |
184 |
> |
guarantees that all of the random forces point inward towards the |
185 |
> |
cluster. |
186 |
|
|
187 |
< |
Novel features: No a priori geometry is defined, No affine transforms, |
188 |
< |
No fictitious particles, No bounding potentials. |
187 |
> |
In the following sections, we extend and generalize the approach of |
188 |
> |
Kohanoff, Caro, and Finnis. The new method, which we are calling the |
189 |
> |
``Langevin Hull'' applies the external pressure, Langevin drag, and |
190 |
> |
random forces on the facets of the {\it hull itself} instead of the |
191 |
> |
atomic sites comprising the vertices of the hull. This allows us to |
192 |
> |
decouple the external pressure contribution from the drag and random |
193 |
> |
force. Section \ref{sec:meth} |
194 |
|
|
195 |
< |
Simulation starts as a collection of atomic locations in 3D (a point |
196 |
< |
cloud). |
195 |
> |
\section{Methodology} |
196 |
> |
\label{sec:meth} |
197 |
|
|
198 |
< |
Delaunay triangulation finds all facets between coplanar neighbors. |
198 |
> |
We have developed a new method which uses a constant pressure and |
199 |
> |
temperature bath. This bath interacts only with the objects that are |
200 |
> |
currently at the edge of the system. Since the edge is determined |
201 |
> |
dynamically as the simulation progresses, no {\it a priori} geometry |
202 |
> |
is defined. The pressure and temperature bath interacts {\it |
203 |
> |
directly} with the atoms on the edge and not with atoms interior to |
204 |
> |
the simulation. This means that there are no affine transforms |
205 |
> |
required. There are also no fictitious particles or bounding |
206 |
> |
potentials used in this approach. |
207 |
|
|
208 |
< |
The Convex Hull is the set of facets that have no concave corners at a |
209 |
< |
vertex. |
208 |
> |
The basics of the method are as follows. The simulation starts as a |
209 |
> |
collection of atomic locations in three dimensions (a point cloud). |
210 |
> |
Delaunay triangulation is used to find all facets between coplanar |
211 |
> |
neighbors. In highly symmetric point clouds, facets can contain many |
212 |
> |
atoms, but in all but the most symmetric of cases one might experience |
213 |
> |
in a molecular dynamics simulation, the facets are simple triangles in |
214 |
> |
3-space that contain exactly three atoms. |
215 |
|
|
216 |
< |
Molecules on the convex hull are dynamic. As they re-enter the |
217 |
< |
cluster, all interactions with the external bath are removed.The |
218 |
< |
external bath applies pressure to the facets of the convex hull in |
219 |
< |
direct proportion to the area of the facet. Thermal coupling depends on |
220 |
< |
the solvent temperature, friction and the size and shape of each |
221 |
< |
facet. |
216 |
> |
The convex hull is the set of facets that have {\it no concave |
217 |
> |
corners} at an atomic site. This eliminates all facets on the |
218 |
> |
interior of the point cloud, leaving only those exposed to the |
219 |
> |
bath. Sites on the convex hull are dynamic. As molecules re-enter the |
220 |
> |
cluster, all interactions between atoms on that molecule and the |
221 |
> |
external bath are removed. |
222 |
|
|
223 |
+ |
For atomic sites in the interior of the point cloud, the equations of |
224 |
+ |
motion are simple Newtonian dynamics, |
225 |
|
\begin{equation} |
226 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U |
226 |
> |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
227 |
> |
\label{eq:Newton} |
228 |
|
\end{equation} |
229 |
< |
|
229 |
> |
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the |
230 |
> |
instantaneous velocity of site $i$ at time $t$, and $U$ is the total |
231 |
> |
potential energy. For atoms on the exterior of the cluster |
232 |
> |
(i.e. those that occupy one of the vertices of the convex hull), the |
233 |
> |
equation of motion is modified with an external force, ${\mathbf |
234 |
> |
F}_i^{\mathrm ext}$, |
235 |
|
\begin{equation} |
236 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext} |
236 |
> |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
237 |
|
\end{equation} |
238 |
|
|
239 |
+ |
The external bath interacts directly with the facets of the convex |
240 |
+ |
hull. Since each vertex (or atom) provides one corner of a triangular |
241 |
+ |
facet, the force on the facets are divided equally to each vertex. |
242 |
+ |
However, each vertex can participate in multiple facets, so the resultant |
243 |
+ |
force is a sum over all facets $f$ containing vertex $i$: |
244 |
|
\begin{equation} |
245 |
|
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\ |
246 |
|
} f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf |
247 |
|
F}_f^{\mathrm ext} |
248 |
|
\end{equation} |
249 |
|
|
250 |
+ |
The external pressure bath applies a force to the facets of the convex |
251 |
+ |
hull in direct proportion to the area of the facet, while the thermal |
252 |
+ |
coupling depends on the solvent temperature, friction and the size and |
253 |
+ |
shape of each facet. The thermal interactions are expressed as a |
254 |
+ |
typical Langevin description of the forces, |
255 |
|
\begin{equation} |
256 |
|
\begin{array}{rclclcl} |
257 |
|
{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\ |
258 |
|
& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t) |
259 |
|
\end{array} |
260 |
|
\end{equation} |
261 |
< |
|
261 |
> |
Here, $P$ is the external pressure, $A_f$ and $\hat{n}_f$ are the area |
262 |
> |
and normal vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is |
263 |
> |
the velocity of the facet, |
264 |
> |
\begin{equation} |
265 |
> |
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
266 |
> |
\end{equation} |
267 |
> |
and $\Xi_f(t)$ is an approximate ($3 \times 3$) hydrodynamic tensor |
268 |
> |
that depends on the geometry and surface area of facet $f$ and the |
269 |
> |
viscosity of the fluid (See Appendix A). The hydrodynamic tensor is |
270 |
> |
related to the fluctuations of the random force, $\mathbf{R}(t)$, by |
271 |
> |
the fluctuation-dissipation theorem, |
272 |
|
\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} |
273 |
|
\left< {\mathbf R}_f(t) \right> & = & 0 \\ |
274 |
|
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\ |
275 |
< |
\Xi_f(t)\delta(t-t^\prime) |
275 |
> |
\Xi_f(t)\delta(t-t^\prime). |
276 |
> |
\label{eq:randomForce} |
277 |
|
\end{eqnarray} |
278 |
|
|
279 |
< |
Implemented in OpenMD.\cite{Meineke2005,openmd} |
279 |
> |
Once the hydrodynamic tensor is known for a given facet (see Appendix |
280 |
> |
A) obtaining a stochastic vector that has the properties in |
281 |
> |
Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a |
282 |
> |
one-time Cholesky decomposition to obtain the square root matrix of |
283 |
> |
the resistance tensor, |
284 |
> |
\begin{equation} |
285 |
> |
\Xi_f = {\bf S} {\bf S}^{T}, |
286 |
> |
\label{eq:Cholesky} |
287 |
> |
\end{equation} |
288 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
289 |
> |
vector with the statistics required for the random force can then be |
290 |
> |
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which |
291 |
> |
has elements chosen from a Gaussian distribution, such that: |
292 |
> |
\begin{equation} |
293 |
> |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
294 |
> |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
295 |
> |
\end{equation} |
296 |
> |
where $\delta t$ is the timestep in use during the simulation. The |
297 |
> |
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
298 |
> |
have the correct properties required by Eq. (\ref{eq:randomForce}). |
299 |
|
|
300 |
+ |
Our treatment of the hydrodynamic tensor must be approximate. $\Xi$ |
301 |
+ |
for a triangular plate would normally be treated as a $6 \times 6$ |
302 |
+ |
tensor that includes translational and rotational drag as well as |
303 |
+ |
translational-rotational coupling. The computation of hydrodynamic |
304 |
+ |
tensors for rigid bodies has been detailed |
305 |
+ |
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun2008} |
306 |
+ |
but the standard approach involving bead approximations would be |
307 |
+ |
prohibitively expensive if it were recomputed at each step in a |
308 |
+ |
molecular dynamics simulation. |
309 |
+ |
|
310 |
+ |
We are utilizing an approximate hydrodynamic tensor obtained by first |
311 |
+ |
constructing the Oseen tensor for the interaction of the centroid of |
312 |
+ |
the facet ($f$) with each of the subfacets $j$, |
313 |
+ |
\begin{equation} |
314 |
+ |
T_{jf}=\frac{A_j}{8\pi\eta R_{jf}}\left(I + |
315 |
+ |
\frac{\mathbf{R}_{jf}\mathbf{R}_{jf}^T}{R_{jf}^2}\right) |
316 |
+ |
\end{equation} |
317 |
+ |
Here, $A_j$ is the area of subfacet $j$ which is a triangle containing |
318 |
+ |
two of the vertices of the facet along with the centroid. |
319 |
+ |
$\mathbf{R}_{jf}$ is the vector between the centroid of facet $f$ and |
320 |
+ |
the centroid of sub-facet $j$, and $I$ is the ($3 \times 3$) identity |
321 |
+ |
matrix. $\eta$ is the viscosity of the external bath. |
322 |
+ |
|
323 |
+ |
\begin{figure} |
324 |
+ |
\includegraphics[width=\linewidth]{hydro} |
325 |
+ |
\caption{The hydrodynamic tensor $\Xi$ for a facet comprising sites $i$, |
326 |
+ |
$j$, and $k$ is constructed using Oseen tensor contributions |
327 |
+ |
between the centoid of the facet $f$ and each of the sub-facets |
328 |
+ |
($i,f,j$), ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets |
329 |
+ |
are located at $1$, $2$, and $3$, and the area of each sub-facet is |
330 |
+ |
easily computed using half the cross product of two of the edges.} |
331 |
+ |
\label{hydro} |
332 |
+ |
\end{figure} |
333 |
+ |
|
334 |
+ |
The Oseen tensors for each of the sub-facets are summed, and the |
335 |
+ |
resulting matrix is inverted to give a $3 \times 3$ hydrodynamic |
336 |
+ |
tensor for translations of the triangular plate, |
337 |
+ |
\begin{equation} |
338 |
+ |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}. |
339 |
+ |
\end{equation} |
340 |
+ |
We have implemented this method by extending the Langevin dynamics |
341 |
+ |
integrator in our group code, OpenMD.\cite{Meineke2005,openmd} There |
342 |
+ |
is a moderate penalty for computing the convex hull at each step in |
343 |
+ |
the molecular dynamics simulation (HOW MUCH?), but the convex hull is |
344 |
+ |
remarkably easy to parallelize on distributed memory machines (see |
345 |
+ |
Appendix B). |
346 |
+ |
|
347 |
|
\section{Tests \& Applications} |
348 |
+ |
\label{sec:tests} |
349 |
|
|
350 |
|
\subsection{Bulk modulus of gold nanoparticles} |
351 |
|
|
452 |
|
|
453 |
|
\section{Appendix A: Hydrodynamic tensor for triangular facets} |
454 |
|
|
278 |
– |
\begin{figure} |
279 |
– |
\includegraphics[width=\linewidth]{hydro} |
280 |
– |
\caption{Hydro} |
281 |
– |
\label{hydro} |
282 |
– |
\end{figure} |
283 |
– |
|
284 |
– |
\begin{equation} |
285 |
– |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1} |
286 |
– |
\end{equation} |
287 |
– |
|
288 |
– |
\begin{equation} |
289 |
– |
T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I + |
290 |
– |
\frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right) |
291 |
– |
\end{equation} |
292 |
– |
|
455 |
|
\section{Appendix B: Computing Convex Hulls on Parallel Computers} |
456 |
|
|
457 |
|
\section{Acknowledgments} |