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. |
84 |
< |
|
85 |
< |
As long as the material in the simulation box is essentially a bulk |
86 |
< |
liquid which has a relatively uniform compressibility, the standard |
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. |
88 |
> |
atomic sites. |
89 |
|
|
90 |
< |
The problem with these approaches becomes apparent when the material |
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 |
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 collisions in the incompressible region. |
99 |
> |
enough to avoid the instabilities in the incompressible region. |
100 |
|
|
101 |
|
\begin{figure} |
102 |
|
\includegraphics[width=\linewidth]{AffineScale2} |
110 |
|
\label{affineScale} |
111 |
|
\end{figure} |
112 |
|
|
113 |
< |
Additionally, one may often wish to simulate explicitly non-periodic |
114 |
< |
systems, and the constraint that a periodic box must be used to |
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 |
< |
Explicitly non-periodic systems |
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 |
< |
Elastic Bag |
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 |
< |
Spherical Boundary approaches |
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 |
+ |
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 |
+ |
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 |
+ |
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 |
+ |
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 |
|
\section{Methodology} |
196 |
+ |
\label{sec:meth} |
197 |
|
|
198 |
|
We have developed a new method which uses a constant pressure and |
199 |
|
temperature bath. This bath interacts only with the objects that are |
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 a ($3 \times 3$) hydrodynamic tensor that depends on |
268 |
< |
the geometry and surface area of facet $f$ and the viscosity of the |
269 |
< |
fluid (See Appendix A). The hydrodynamic tensor is related to the |
270 |
< |
fluctuations of the random force, $\mathbf{R}(t)$, by the |
271 |
< |
fluctuation-dissipation theorem, |
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} |
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\ |
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} |
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 |
|
|
337 |
– |
\begin{figure} |
338 |
– |
\includegraphics[width=\linewidth]{hydro} |
339 |
– |
\caption{Hydro} |
340 |
– |
\label{hydro} |
341 |
– |
\end{figure} |
342 |
– |
|
343 |
– |
\begin{equation} |
344 |
– |
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1} |
345 |
– |
\end{equation} |
346 |
– |
|
347 |
– |
\begin{equation} |
348 |
– |
T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I + |
349 |
– |
\frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right) |
350 |
– |
\end{equation} |
351 |
– |
|
455 |
|
\section{Appendix B: Computing Convex Hulls on Parallel Computers} |
456 |
|
|
457 |
|
\section{Acknowledgments} |