| 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} |
