| 24 |
|
|
| 25 |
|
\title{The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems} |
| 26 |
|
|
| 27 |
< |
\author{Charles F. Varedeman II, Kelsey Stocker, and J. Daniel |
| 27 |
> |
\author{Charles F. Vardeman II, Kelsey M. Stocker, and J. Daniel |
| 28 |
|
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ |
| 29 |
|
Department of Chemistry and Biochemistry,\\ |
| 30 |
|
University of Notre Dame\\ |
| 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 |
| 65 |
|
|
| 66 |
|
\section{Introduction} |
| 67 |
|
|
| 68 |
< |
Affine transform methods |
| 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 |
| 77 |
> |
particles). The most common constant pressure methods, including the |
| 78 |
> |
Melchionna modification\cite{Melchionna1993} to the |
| 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 |
< |
\begin{figure} |
| 86 |
< |
\includegraphics[width=\linewidth]{AffineScale} |
| 87 |
< |
\caption{Affine Scale} |
| 88 |
< |
\label{affineScale} |
| 89 |
< |
\end{figure} |
| 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 Scale2} |
| 105 |
< |
\label{affineScale2} |
| 104 |
> |
\caption{Affine Scaling constant pressure methods use box-length |
| 105 |
> |
scaling to adjust the volume to adjust to under- or over-pressure |
| 106 |
> |
conditions. In a system with a uniform compressibility (e.g. bulk |
| 107 |
> |
fluids) these methods can work well. In systems containing |
| 108 |
> |
heterogeneous mixtures, the affine scaling moves required to adjust |
| 109 |
> |
the pressure in the high-compressibility regions can cause molecules |
| 110 |
> |
in low compressibility regions to collide.} |
| 111 |
> |
\label{affineScale} |
| 112 |
|
\end{figure} |
| 113 |
|
|
| 114 |
< |
Heterogeneous mixtures of materials with different compressibilities? |
| 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 |
|
|
| 117 |
|
Explicitly non-periodic systems |
| 118 |
|
|
| 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). |
| 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 |
< |
Delaunay triangulation finds all facets between coplanar neighbors. |
| 151 |
< |
|
| 105 |
< |
The Convex Hull is the set of facets that have no concave corners at a |
| 106 |
< |
vertex. |
| 107 |
< |
|
| 108 |
< |
Molecules on the convex hull are dynamic. As they re-enter the |
| 109 |
< |
cluster, all interactions with the external bath are removed.The |
| 110 |
< |
external bath applies pressure to the facets of the convex hull in |
| 111 |
< |
direct proportion to the area of the facet.Thermal coupling depends on |
| 112 |
< |
the solvent temperature, friction and the size and shape of each |
| 113 |
< |
facet. |
| 114 |
< |
|
| 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} |
| 137 |
– |
A_f & = & \text{area of facet}\ f \\ |
| 138 |
– |
\hat{n}_f & = & \text{facet normal} \\ |
| 139 |
– |
P & = & \text{external pressure} |
| 140 |
– |
\end{eqnarray} |
| 141 |
– |
|
| 142 |
– |
\begin{eqnarray} |
| 143 |
– |
{\mathbf v}_f(t) & = & \text{velocity of facet} \\ |
| 144 |
– |
& = & \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i \\ |
| 145 |
– |
\Xi_f(t) & = & \text{is a hydrodynamic tensor that depends} \\ |
| 146 |
– |
& & \text{on the geometry and surface area of} \\ |
| 147 |
– |
& & \text{facet} \ f\ \text{and the viscosity of the fluid.} |
| 148 |
– |
\end{eqnarray} |
| 149 |
– |
|
| 150 |
– |
\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{Meineke:2005gd,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} |
| 258 |
|
|
| 259 |
|
\subsection{Compressibility of SPC/E water clusters} |
| 260 |
|
|
| 261 |
+ |
Both NVT \cite{Glattli2002} and NPT \cite{Motakabbir1990, Pi2009} molecular dynamics simulations of SPC/E water have yielded values for the isothermal compressibility of water that agree well with experiment \cite{Fine1973}. The results of three different methods for computing the isothermal compressibility from Langevin Hull simulations for pressures between 1 and 6500 atm are shown in Fig. 5 along with compressibility values obtained from both other SPC/E simulations and experiment. Compressibility values from all references are for applied pressures within the range 1 - 1000 atm. |
| 262 |
+ |
|
| 263 |
|
\begin{figure} |
| 264 |
< |
\includegraphics[width=\linewidth]{g_r_theta} |
| 265 |
< |
\caption{Definition of coordinates} |
| 266 |
< |
\label{coords} |
| 264 |
> |
\includegraphics[width=\linewidth]{new_isothermal} |
| 265 |
> |
\caption{Compressibility of SPC/E water} |
| 266 |
> |
\label{compWater} |
| 267 |
|
\end{figure} |
| 268 |
|
|
| 269 |
+ |
We initially used the classic compressibility formula |
| 270 |
+ |
|
| 271 |
|
\begin{equation} |
| 272 |
+ |
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T} |
| 273 |
+ |
\end{equation} |
| 274 |
+ |
|
| 275 |
+ |
to calculate the the isothermal compressibility at each target pressure. These calculations yielded compressibility values that were dramatically higher than both previous simulations and experiment. The particular compressibility expression used requires the calculation of both a volume and pressure differential, thereby stipulating that the data from at least two simulations at different pressures must be used to calculate the isothermal compressibility at one pressure. |
| 276 |
+ |
|
| 277 |
+ |
Per the fluctuation dissipation theorem \cite{Debenedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure |
| 278 |
+ |
|
| 279 |
+ |
\begin{equation} |
| 280 |
+ |
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T} |
| 281 |
+ |
\end{equation} |
| 282 |
+ |
|
| 283 |
+ |
Thus, the compressibility of each simulation run can be calculated entirely independently from all other trajectories. However, the resulting compressibilities were still as much as an order of magnitude larger than the reference values. The effect was particularly pronounced at the low end of the pressure range. At ambient temperature and low pressures, there exists an equilibrium between vapor and liquid phases. Vapor molecules are naturally more diffuse around the exterior of the cluster, causing artificially large cluster volumes. Any compressibility calculation that relies on the hull volume will suffer these effects. |
| 284 |
+ |
|
| 285 |
+ |
In order to calculate the isothermal compressibility without being hindered by hull volume issues, we adapted the classic compressibility formula so that the compressibility could be calculated using information about the local density instead of the volume of the convex hull. We calculated the $g_{OO}(r)$ for a 1 nanosecond simulation of a cluster of 1372 SPC/E water molecules and spherically integrated the function over the bounds 0 to $r'$. In all cases, the value of $r'$ was 17.26216 $\AA$. The value of the total integral between these bounds is essentially the number (N) of molecules within volume $\frac{4}{3}\pi r'^{3}$ at a given pressure. To yield an actual molecule count, N must be scaled by an ideal density. However, even in the absence of an ideal density, we can use the relationship $\rho = \frac{N}{V}$ to rewrite the isothermal compressibility formula as |
| 286 |
+ |
|
| 287 |
+ |
\begin{equation} |
| 288 |
+ |
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T} |
| 289 |
+ |
\end{equation} |
| 290 |
+ |
|
| 291 |
+ |
Isothermal compressibility values calculated using this modified expression are in good agreement with the reference values throughout the 1 - 1000 atm pressure regime. Regardless of the difficulty in obtaining accurate hull volumes at low temperature and pressures, the Langevin Hull NPT method provides reasonable isothermal compressibility values for water through a large range of pressures. |
| 292 |
+ |
|
| 293 |
+ |
\subsection{Molecular orientation distribution at cluster boundary} |
| 294 |
+ |
|
| 295 |
+ |
In order for non-periodic boundary conditions to be widely applicable, they must be constructed in such a way that they allow a finite, usually small, simulated system to replicate the properties of an infinite bulk system. Naturally, this requirement has spawned many methods for inserting boundaries into simulated systems [REF... ?]. Of particular interest to our characterization of the Langevin Hull is the orientation of water molecules included in the geometric hull. Ideally, all molecules in the cluster will have the same orientational distribution as bulk water. |
| 296 |
+ |
|
| 297 |
+ |
The orientation of molecules at the edges of a simulated cluster has long been a concern when performing simulations of explicitly non-periodic systems. Early work led to the surface constrained soft sphere dipole model (SCSSD) \cite{Warshel1978} in which the surface molecules are fixed in a random orientation representative of the bulk solvent structural properties. Belch, et al \cite{Belch1985} simulated clusters of TIPS2 water surrounded by a hydrophobic bounding potential. The spherical hydrophobic boundary induced dangling hydrogen bonds at the surface that propagated deep into the cluster, affecting 70\% of the 100 molecules in the simulation. This result echoes an earlier study which showed that an extended planar hydrophobic surface caused orientational preference at the surface which extended 7 \r{A} into the liquid simulation cell \cite{Lee1984}. The surface constrained all-atom solvent (SCAAS) model \cite{King1989} improved upon its SCSSD predecessor. The SCAAS model utilizes a polarization constraint which is applied to the surface molecules to maintain bulk-like structure at the cluster surface. A radial constraint is used to maintain the desired bulk density of the liquid. Both constraint forces are applied only to a pre-determined number of the outermost molecules. |
| 298 |
+ |
|
| 299 |
+ |
In contrast, the Langevin Hull does not require that the orientation of molecules be fixed, nor does it utilize an explicitly hydrophobic boundary, orientational constraint or radial constraint. The number and identity of the molecules included on the convex hull are dynamic properties, thus avoiding the formation of an artificial solvent boundary layer. The hope is that the water molecules on the surface of the cluster, if left to their own devices in the absence of orientational and radial constraints, will maintain a bulk-like orientational distribution. |
| 300 |
+ |
|
| 301 |
+ |
To determine the extent of these effects demonstrated by the Langevin Hull, we examined the orientations exhibited by SPC/E water in a cluster of 1372 molecules at 300 K and at pressures ranging from 1 - 1000 atm. |
| 302 |
+ |
|
| 303 |
+ |
The orientation of a water molecule is described by |
| 304 |
+ |
|
| 305 |
+ |
\begin{equation} |
| 306 |
|
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|} |
| 307 |
|
\end{equation} |
| 308 |
|
|
| 309 |
+ |
where $\vec{r}_{i}$ is the vector between molecule {\it i}'s center of mass and the cluster center of mass and $\vec{\mu}_{i}$ is the vector bisecting the H-O-H angle of molecule {\it i}. |
| 310 |
+ |
|
| 311 |
|
\begin{figure} |
| 312 |
+ |
\includegraphics[width=\linewidth]{g_r_theta} |
| 313 |
+ |
\caption{Definition of coordinates} |
| 314 |
+ |
\label{coords} |
| 315 |
+ |
\end{figure} |
| 316 |
+ |
|
| 317 |
+ |
Fig. 7 shows the probability of each value of $\cos{\theta}$ for molecules in the interior of the cluster (squares) and for molecules included in the convex hull (circles). |
| 318 |
+ |
|
| 319 |
+ |
\begin{figure} |
| 320 |
|
\includegraphics[width=\linewidth]{pAngle} |
| 321 |
|
\caption{SPC/E water clusters: only minor dewetting at the boundary} |
| 322 |
|
\label{pAngle} |
| 323 |
|
\end{figure} |
| 324 |
|
|
| 325 |
< |
\begin{figure} |
| 206 |
< |
\includegraphics[width=\linewidth]{isothermal} |
| 207 |
< |
\caption{Compressibility of SPC/E water} |
| 208 |
< |
\label{compWater} |
| 209 |
< |
\end{figure} |
| 325 |
> |
As expected, interior molecules (those not included in the convex hull) maintain a bulk-like structure with a uniform distribution of orientations. Molecules included in the convex hull show a slight preference for values of $\cos{\theta} < 0.$ These values correspond to molecules with a hydrogen directed toward the exterior of the cluster, forming a dangling hydrogen bond. |
| 326 |
|
|
| 327 |
+ |
In the absence of an electrostatic contribution from the exterior bath, the orientational distribution of water molecules included in the Langevin Hull will slightly resemble the distribution at a neat water liquid/vapor interface. Previous molecular dynamics simulations of SPC/E water \cite{Taylor1996} have shown that molecules at the liquid/vapor interface favor an orientation where one hydrogen protrudes from the liquid phase. This behavior is demonstrated by experiments \cite{Du1994} \cite{Scatena2001} showing that approximately one-quarter of water molecules at the liquid/vapor interface form dangling hydrogen bonds. The negligible preference shown in these cluster simulations could be removed through the introduction of an implicit solvent model, which would provide the missing electrostatic interactions between the cluster molecules and the surrounding temperature/pressure bath. |
| 328 |
+ |
|
| 329 |
+ |
The orientational preference exhibited by hull molecules is significantly weaker than the preference caused by an explicit hydrophobic bounding potential. Additionally, the Langevin Hull does not require that the orientation of any molecules be fixed in order to maintain bulk-like structure, even at the cluster surface. |
| 330 |
+ |
|
| 331 |
+ |
|
| 332 |
|
\subsection{Heterogeneous nanoparticle / water mixtures} |
| 333 |
|
|
| 334 |
|
|
