| 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. |
| 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 |
| 87 |
+ |
approach provides an excellent way of adjusting the volume of the |
| 88 |
+ |
system and applying pressure directly via the interactions between |
| 89 |
+ |
atomic sites. |
| 90 |
+ |
|
| 91 |
+ |
The problem with these approaches becomes apparent when the material |
| 92 |
+ |
being simulated is an inhomogeneous mixture in which portions of the |
| 93 |
+ |
simulation box are incompressible relative to other portions. |
| 94 |
+ |
Examples include simulations of metallic nanoparticles in liquid |
| 95 |
+ |
environments, proteins at interfaces, as well as other multi-phase or |
| 96 |
+ |
interfacial environments. In these cases, the affine transform of |
| 97 |
+ |
atomic coordinates will either cause numerical instability when the |
| 98 |
+ |
sites in the incompressible medium collide with each other, or lead to |
| 99 |
+ |
inefficient sampling of system volumes if the barostat is set slow |
| 100 |
+ |
enough to avoid collisions in the incompressible region. |
| 101 |
+ |
|
| 102 |
|
\begin{figure} |
| 103 |
|
\includegraphics[width=\linewidth]{AffineScale2} |
| 104 |
|
\caption{Affine Scaling constant pressure methods use box-length |
| 111 |
|
\label{affineScale} |
| 112 |
|
\end{figure} |
| 113 |
|
|
| 114 |
+ |
Additionally, one may often wish to simulate explicitly non-periodic |
| 115 |
+ |
systems, and the constraint that a periodic box must be used to |
| 116 |
|
|
| 96 |
– |
Heterogeneous mixtures of materials with different compressibilities? |
| 97 |
– |
|
| 117 |
|
Explicitly non-periodic systems |
| 118 |
|
|
| 119 |
|
Elastic Bag |
| 122 |
|
|
| 123 |
|
\section{Methodology} |
| 124 |
|
|
| 125 |
< |
A new method which uses a constant pressure and temperature bath that |
| 126 |
< |
interacts with the objects that are currently at the edge of the |
| 127 |
< |
system. |
| 125 |
> |
We have developed a new method which uses a constant pressure and |
| 126 |
> |
temperature bath. This bath interacts only with the objects that are |
| 127 |
> |
currently at the edge of the system. Since the edge is determined |
| 128 |
> |
dynamically as the simulation progresses, no {\it a priori} geometry |
| 129 |
> |
is defined. The pressure and temperature bath interacts {\it |
| 130 |
> |
directly} with the atoms on the edge and not with atoms interior to |
| 131 |
> |
the simulation. This means that there are no affine transforms |
| 132 |
> |
required. There are also no fictitious particles or bounding |
| 133 |
> |
potentials used in this approach. |
| 134 |
|
|
| 135 |
< |
Novel features: No a priori geometry is defined, No affine transforms, |
| 136 |
< |
No fictitious particles, No bounding potentials. |
| 135 |
> |
The basics of the method are as follows. The simulation starts as a |
| 136 |
> |
collection of atomic locations in three dimensions (a point cloud). |
| 137 |
> |
Delaunay triangulation is used to find all facets between coplanar |
| 138 |
> |
neighbors. In highly symmetric point clouds, facets can contain many |
| 139 |
> |
atoms, but in all but the most symmetric of cases one might experience |
| 140 |
> |
in a molecular dynamics simulation, the facets are simple triangles in |
| 141 |
> |
3-space that contain exactly three atoms. |
| 142 |
|
|
| 143 |
< |
Simulation starts as a collection of atomic locations in 3D (a point |
| 144 |
< |
cloud). |
| 145 |
< |
|
| 146 |
< |
Delaunay triangulation finds all facets between coplanar neighbors. |
| 147 |
< |
|
| 148 |
< |
The Convex Hull is the set of facets that have no concave corners at a |
| 119 |
< |
vertex. |
| 120 |
< |
|
| 121 |
< |
Molecules on the convex hull are dynamic. As they re-enter the |
| 122 |
< |
cluster, all interactions with the external bath are removed.The |
| 123 |
< |
external bath applies pressure to the facets of the convex hull in |
| 124 |
< |
direct proportion to the area of the facet. Thermal coupling depends on |
| 125 |
< |
the solvent temperature, friction and the size and shape of each |
| 126 |
< |
facet. |
| 143 |
> |
The convex hull is the set of facets that have {\it no concave |
| 144 |
> |
corners} at an atomic site. This eliminates all facets on the |
| 145 |
> |
interior of the point cloud, leaving only those exposed to the |
| 146 |
> |
bath. Sites on the convex hull are dynamic. As molecules re-enter the |
| 147 |
> |
cluster, all interactions between atoms on that molecule and the |
| 148 |
> |
external bath are removed. |
| 149 |
|
|
| 150 |
+ |
For atomic sites in the interior of the point cloud, the equations of |
| 151 |
+ |
motion are simple Newtonian dynamics, |
| 152 |
|
\begin{equation} |
| 153 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U |
| 153 |
> |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U, |
| 154 |
> |
\label{eq:Newton} |
| 155 |
|
\end{equation} |
| 156 |
< |
|
| 156 |
> |
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the |
| 157 |
> |
instantaneous velocity of site $i$ at time $t$, and $U$ is the total |
| 158 |
> |
potential energy. For atoms on the exterior of the cluster |
| 159 |
> |
(i.e. those that occupy one of the vertices of the convex hull), the |
| 160 |
> |
equation of motion is modified with an external force, ${\mathbf |
| 161 |
> |
F}_i^{\mathrm ext}$, |
| 162 |
|
\begin{equation} |
| 163 |
< |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext} |
| 163 |
> |
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}. |
| 164 |
|
\end{equation} |
| 165 |
|
|
| 166 |
+ |
The external bath interacts directly with the facets of the convex |
| 167 |
+ |
hull. Since each vertex (or atom) provides one corner of a triangular |
| 168 |
+ |
facet, the force on the facets are divided equally to each vertex. |
| 169 |
+ |
However, each vertex can participate in multiple facets, so the resultant |
| 170 |
+ |
force is a sum over all facets $f$ containing vertex $i$: |
| 171 |
|
\begin{equation} |
| 172 |
|
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\ |
| 173 |
|
} f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf |
| 174 |
|
F}_f^{\mathrm ext} |
| 175 |
|
\end{equation} |
| 176 |
|
|
| 177 |
+ |
The external pressure bath applies a force to the facets of the convex |
| 178 |
+ |
hull in direct proportion to the area of the facet, while the thermal |
| 179 |
+ |
coupling depends on the solvent temperature, friction and the size and |
| 180 |
+ |
shape of each facet. The thermal interactions are expressed as a |
| 181 |
+ |
typical Langevin description of the forces, |
| 182 |
|
\begin{equation} |
| 183 |
|
\begin{array}{rclclcl} |
| 184 |
|
{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\ |
| 185 |
|
& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t) |
| 186 |
|
\end{array} |
| 187 |
|
\end{equation} |
| 188 |
< |
|
| 188 |
> |
Here, $P$ is the external pressure, $A_f$ and $\hat{n}_f$ are the area |
| 189 |
> |
and normal vectors for facet $f$, respectively. ${\mathbf v}_f(t)$ is |
| 190 |
> |
the velocity of the facet, |
| 191 |
> |
\begin{equation} |
| 192 |
> |
{\mathbf v}_f(t) = \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i, |
| 193 |
> |
\end{equation} |
| 194 |
> |
and $\Xi_f(t)$ is a ($3 \times 3$) hydrodynamic tensor that depends on |
| 195 |
> |
the geometry and surface area of facet $f$ and the viscosity of the |
| 196 |
> |
fluid (See Appendix A). The hydrodynamic tensor is related to the |
| 197 |
> |
fluctuations of the random force, $\mathbf{R}(t)$, by the |
| 198 |
> |
fluctuation-dissipation theorem, |
| 199 |
|
\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} |
| 200 |
|
\left< {\mathbf R}_f(t) \right> & = & 0 \\ |
| 201 |
|
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\ |
| 202 |
< |
\Xi_f(t)\delta(t-t^\prime) |
| 202 |
> |
\Xi_f(t)\delta(t-t^\prime). |
| 203 |
> |
\label{eq:randomForce} |
| 204 |
|
\end{eqnarray} |
| 205 |
|
|
| 206 |
< |
Implemented in OpenMD.\cite{Meineke2005,openmd} |
| 206 |
> |
Once the hydrodynamic tensor is known for a given facet (see Appendix |
| 207 |
> |
A) obtaining a stochastic vector that has the properties in |
| 208 |
> |
Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a |
| 209 |
> |
one-time Cholesky decomposition to obtain the square root matrix of |
| 210 |
> |
the resistance tensor, |
| 211 |
> |
\begin{equation} |
| 212 |
> |
\Xi_f = {\bf S} {\bf S}^{T}, |
| 213 |
> |
\label{eq:Cholesky} |
| 214 |
> |
\end{equation} |
| 215 |
> |
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A |
| 216 |
> |
vector with the statistics required for the random force can then be |
| 217 |
> |
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which |
| 218 |
> |
has elements chosen from a Gaussian distribution, such that: |
| 219 |
> |
\begin{equation} |
| 220 |
> |
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot |
| 221 |
> |
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, |
| 222 |
> |
\end{equation} |
| 223 |
> |
where $\delta t$ is the timestep in use during the simulation. The |
| 224 |
> |
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to |
| 225 |
> |
have the correct properties required by Eq. (\ref{eq:randomForce}). |
| 226 |
|
|
| 227 |
+ |
We have implemented this method by extending the Langevin dynamics |
| 228 |
+ |
integrator in our group code, OpenMD.\cite{Meineke2005,openmd} |
| 229 |
+ |
|
| 230 |
|
\section{Tests \& Applications} |
| 231 |
|
|
| 232 |
|
\subsection{Bulk modulus of gold nanoparticles} |
