trunk/langevinHull/langevinHull.tex

\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
\usepackage{graphicx}
\usepackage{multirow}
\usepackage[square, comma, sort&compress]{natbib}
\usepackage{url}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000

% double space list of tables and figures
%\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

\bibpunct{[}{]}{,}{s}{}{;}
\bibliographystyle{aip}

\begin{document}

\title{The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems}

\author{Charles F. Vardeman II, Kelsey M. Stocker, and J. Daniel
Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
Department of Chemistry and Biochemistry,\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle 

\begin{doublespace}

\begin{abstract}
  We have developed a new isobaric-isothermal (NPT) algorithm which
  applies an external pressure to the facets comprising the convex
  hull surrounding the objects in the system. Additionally, a Langevin
  thermostat is applied to facets of the hull to mimic contact with an
  external heat bath. This new method, the ``Langevin Hull'',
  performs better than traditional affine transform methods for
  systems containing heterogeneous mixtures of materials with
  different compressibilities. It does not suffer from the edge
  effects of boundary potential methods, and allows realistic
  treatment of both external pressure and thermal conductivity to an
  implicit solvents.  We apply this method to several different
  systems including bare nano-particles, nano-particles in explicit
  solvent, as well as clusters of liquid water and ice. The predicted
  mechanical and thermal properties of these systems are in good
  agreement with experimental data.
\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\section{Introduction}

The most common molecular dynamics methods for sampling configurations
of an isobaric-isothermal (NPT) ensemble attempt to maintain a target
pressure in a simulation by coupling the volume of the system to an
extra degree of freedom, the {\it barostat}.  These methods require
periodic boundary conditions, because when the instantaneous pressure
in the system differs from the target pressure, the volume is
typically reduced or expanded using {\it affine transforms} of the
system geometry. An affine transform scales both the box lengths as
well as the scaled particle positions (but not the sizes of the
particles). The most common constant pressure methods, including the
Melchionna modification\cite{melchionna93} to the
Nos\'e-Hoover-Andersen equations of motion, the Berendsen pressure
bath, and the Langevin Piston, all utilize coordinate transformation
to adjust the box volume.

\begin{figure}
\includegraphics[width=\linewidth]{AffineScale2}
\caption{Affine Scaling constant pressure methods use box-length
  scaling to adjust the volume to adjust to under- or over-pressure
  conditions. In a system with a uniform compressibility (e.g. bulk
  fluids) these methods can work well.  In systems containing
  heterogeneous mixtures, the affine scaling moves required to adjust
  the pressure in the high-compressibility regions can cause molecules
  in low compressibility regions to collide.}
\label{affineScale}
\end{figure}


Heterogeneous mixtures of materials with different compressibilities?

Explicitly non-periodic systems

Elastic Bag 

Spherical Boundary approaches

\section{Methodology}

A new method which uses a constant pressure and temperature bath that
interacts with the objects that are currently at the edge of the
system.

Novel features: No a priori geometry is defined, No affine transforms,
No fictitious particles, No bounding potentials.

Simulation starts as a collection of atomic locations in 3D (a point
cloud).

Delaunay triangulation finds all facets between coplanar neighbors.

The Convex Hull is the set of facets that have no concave corners at a
vertex.

Molecules on the convex hull are dynamic. As they re-enter the
cluster, all interactions with the external bath are removed.The
external bath applies pressure to the facets of the convex hull in
direct proportion to the area of the facet. Thermal coupling depends on
the solvent temperature, friction and the size and shape of each
facet. 

\begin{equation}
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U
\end{equation}

\begin{equation}
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}
\end{equation}

\begin{equation}
{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\
    } f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\  {\mathbf
  F}_f^{\mathrm ext}
\end{equation}

\begin{equation}
\begin{array}{rclclcl}
{\mathbf F}_f^{\text{ext}} & = &  \text{external pressure} & + & \text{drag force} & + & \text{random force} \\
& = &  -\hat{n}_f P A_f  & - & \Xi_f(t) {\mathbf v}_f(t)  & + & {\mathbf R}_f(t)
\end{array}
\end{equation}

\begin{eqnarray}
A_f & = & \text{area of facet}\ f \\
\hat{n}_f & = & \text{facet normal} \\
P & = & \text{external pressure}
\end{eqnarray}

\begin{eqnarray}
{\mathbf v}_f(t) & = & \text{velocity of facet} \\
 & = & \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i \\
\Xi_f(t) & = & \text{is a hydrodynamic tensor that depends} \\ 
& & \text{on the geometry and surface area of} \\
& & \text{facet} \ f\ \text{and the viscosity of the fluid.}
\end{eqnarray}

\begin{eqnarray}
\left< {\mathbf R}_f(t) \right> & = & 0 \\
\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\
\Xi_f(t)\delta(t-t^\prime)
\end{eqnarray}

Implemented in OpenMD.\cite{Meineke:2005gd,openmd}

\section{Tests \& Applications}

\subsection{Bulk modulus of gold nanoparticles}

\begin{figure}
\includegraphics[width=\linewidth]{pressure_tb}
\caption{Pressure response is rapid (18 \AA gold nanoparticle), target
pressure = 4 GPa}
\label{pressureResponse}
\end{figure}

\begin{figure}
\includegraphics[width=\linewidth]{temperature_tb}
\caption{Temperature equilibration depends on surface area and bath
  viscosity.  Target Temperature = 300K}
\label{temperatureResponse}
\end{figure}

\begin{equation}
\kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial
    P}\right)
\end{equation}

\begin{figure}
\includegraphics[width=\linewidth]{compress_tb}
\caption{Isothermal Compressibility (18 \AA gold nanoparticle)}
\label{temperatureResponse}
\end{figure}

\subsection{Compressibility of SPC/E water clusters}

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. 

\begin{figure}
\includegraphics[width=\linewidth]{new_isothermal}
\caption{Compressibility of SPC/E water}
\label{compWater}
\end{figure}

We initially used the classic compressibility formula 

\begin{equation} 
\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T}
\end{equation}

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.

Per the fluctuation dissipation theorem \cite{Debendedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure 

\begin{equation}
\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T}
\end{equation}

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.

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

\begin{equation}
\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T}
\end{equation}

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.

\subsection{Molecular orientation distribution at cluster boundary}

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. 

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.

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.

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.

The orientation of a water molecule is described by

\begin{equation} 
\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{|\vec{r}_i||\vec{\mu}_i|}
\end{equation}

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

\begin{figure}
\includegraphics[width=\linewidth]{g_r_theta}
\caption{Definition of coordinates}
\label{coords}
\end{figure}

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). 

\begin{figure}
\includegraphics[width=\linewidth]{pAngle}
\caption{SPC/E water clusters: only minor dewetting at the boundary}
\label{pAngle}
\end{figure}

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.

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.

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. 


\subsection{Heterogeneous nanoparticle / water mixtures}


\section{Appendix A: Hydrodynamic tensor for triangular facets}

\begin{figure}
\includegraphics[width=\linewidth]{hydro}
\caption{Hydro}
\label{hydro}
\end{figure}

\begin{equation}
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}
\end{equation}

\begin{equation}
T_{if}=\frac{A_i}{8\pi\eta R_{if}}\left(I +
  \frac{\mathbf{R}_{if}\mathbf{R}_{if}^T}{R_{if}^2}\right)
\end{equation}

\section{Appendix B: Computing Convex Hulls on Parallel Computers} 

\section{Acknowledgments}
Support for this project was provided by the
National Science Foundation under grant CHE-0848243. Computational
time was provided by the Center for Research Computing (CRC) at the
University of Notre Dame.  

\newpage

\bibliography{langevinHull}

\end{doublespace}
\end{document}
Revision:	3649
Committed:	Thu Sep 23 21:42:37 2010 UTC (13 years, 9 months ago) by kstocke1
Content type:	application/x-tex
File size:	16593 byte(s)
Log Message:	M langevinHull.bib M langevinHull.tex
#	User	Rev	Content
1	gezelter	3640	\documentclass[11pt]{article}
2			\usepackage{amsmath}
3			\usepackage{amssymb}
4			\usepackage{setspace}
5			\usepackage{endfloat}
6			\usepackage{caption}
7			\usepackage{graphicx}
8			\usepackage{multirow}
9			\usepackage[square, comma, sort&compress]{natbib}
10			\usepackage{url}
11			\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
12			\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
13			9.0in \textwidth 6.5in \brokenpenalty=10000
14
15			% double space list of tables and figures
16			%\AtBeginDelayedFloats{\renewcomand{\baselinestretch}{1.66}}
17			\setlength{\abovecaptionskip}{20 pt}
18			\setlength{\belowcaptionskip}{30 pt}
19
20			\bibpunct{[}{]}{,}{s}{}{;}
21			\bibliographystyle{aip}
22
23			\begin{document}
24
25			\title{The Langevin Hull: Constant pressure and temperature dynamics for non-periodic systems}
26
27	kstocke1	3644	\author{Charles F. Vardeman II, Kelsey M. Stocker, and J. Daniel
28	gezelter	3640	Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\
29			Department of Chemistry and Biochemistry,\\
30			University of Notre Dame\\
31			Notre Dame, Indiana 46556}
32
33			\date{\today}
34
35			\maketitle
36
37			\begin{doublespace}
38
39			\begin{abstract}
40			We have developed a new isobaric-isothermal (NPT) algorithm which
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.
55			\end{abstract}
56
57			\newpage
58
59			%\narrowtext
60
61			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62			% BODY OF TEXT
63			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64
65
66			\section{Introduction}
67
68	gezelter	3641	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{melchionna93} 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.
82	gezelter	3640
83			\begin{figure}
84	gezelter	3641	\includegraphics[width=\linewidth]{AffineScale2}
85			\caption{Affine Scaling constant pressure methods use box-length
86			scaling to adjust the volume to adjust to under- or over-pressure
87			conditions. In a system with a uniform compressibility (e.g. bulk
88			fluids) these methods can work well. In systems containing
89			heterogeneous mixtures, the affine scaling moves required to adjust
90			the pressure in the high-compressibility regions can cause molecules
91			in low compressibility regions to collide.}
92	gezelter	3640	\label{affineScale}
93			\end{figure}
94
95
96			Heterogeneous mixtures of materials with different compressibilities?
97
98			Explicitly non-periodic systems
99
100			Elastic Bag
101
102			Spherical Boundary approaches
103
104			\section{Methodology}
105
106			A new method which uses a constant pressure and temperature bath that
107			interacts with the objects that are currently at the edge of the
108			system.
109
110			Novel features: No a priori geometry is defined, No affine transforms,
111			No fictitious particles, No bounding potentials.
112
113			Simulation starts as a collection of atomic locations in 3D (a point
114			cloud).
115
116			Delaunay triangulation finds all facets between coplanar neighbors.
117
118			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	kstocke1	3649	direct proportion to the area of the facet. Thermal coupling depends on
125	gezelter	3640	the solvent temperature, friction and the size and shape of each
126			facet.
127
128			\begin{equation}
129			m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U
130			\end{equation}
131
132			\begin{equation}
133			m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}
134			\end{equation}
135
136			\begin{equation}
137			{\mathbf F}_{i}^{\mathrm ext} = \sum_{\begin{array}{c}\mathrm{facets\
138			} f \\ \mathrm{containing\ } i\end{array}} \frac{1}{3}\ {\mathbf
139			F}_f^{\mathrm ext}
140			\end{equation}
141
142			\begin{equation}
143			\begin{array}{rclclcl}
144			{\mathbf F}_f^{\text{ext}} & = & \text{external pressure} & + & \text{drag force} & + & \text{random force} \\
145			& = & -\hat{n}_f P A_f & - & \Xi_f(t) {\mathbf v}_f(t) & + & {\mathbf R}_f(t)
146			\end{array}
147			\end{equation}
148
149			\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}
164			\left< {\mathbf R}_f(t) \right> & = & 0 \\
165			\left<{\mathbf R}_f(t) {\mathbf R}_f^T(t^\prime)\right> & = & 2 k_B T\
166			\Xi_f(t)\delta(t-t^\prime)
167			\end{eqnarray}
168
169			Implemented in OpenMD.\cite{Meineke:2005gd,openmd}
170
171			\section{Tests \& Applications}
172
173			\subsection{Bulk modulus of gold nanoparticles}
174
175			\begin{figure}
176			\includegraphics[width=\linewidth]{pressure_tb}
177			\caption{Pressure response is rapid (18 \AA gold nanoparticle), target
178			pressure = 4 GPa}
179			\label{pressureResponse}
180			\end{figure}
181
182			\begin{figure}
183			\includegraphics[width=\linewidth]{temperature_tb}
184			\caption{Temperature equilibration depends on surface area and bath
185			viscosity. Target Temperature = 300K}
186			\label{temperatureResponse}
187			\end{figure}
188
189			\begin{equation}
190			\kappa_T=-\frac{1}{V_{\mathrm{eq}}}\left(\frac{\partial V}{\partial
191			P}\right)
192			\end{equation}
193
194			\begin{figure}
195			\includegraphics[width=\linewidth]{compress_tb}
196			\caption{Isothermal Compressibility (18 \AA gold nanoparticle)}
197			\label{temperatureResponse}
198			\end{figure}
199
200			\subsection{Compressibility of SPC/E water clusters}
201
202	kstocke1	3649	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.
203
204	gezelter	3640	\begin{figure}
205	kstocke1	3649	\includegraphics[width=\linewidth]{new_isothermal}
206			\caption{Compressibility of SPC/E water}
207			\label{compWater}
208	gezelter	3640	\end{figure}
209
210	kstocke1	3649	We initially used the classic compressibility formula
211
212	gezelter	3640	\begin{equation}
213	kstocke1	3649	\kappa_{T} = -\frac{1}{V} \left ( \frac{\partial V}{\partial P} \right )_{T}
214			\end{equation}
215
216			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.
217
218			Per the fluctuation dissipation theorem \cite{Debendedetti1986}, the hull volume fluctuation in any given simulation can be used to calculated the isothermal compressibility at that particular pressure
219
220			\begin{equation}
221			\kappa_{T} = \frac{\left \langle V^{2} \right \rangle - \left \langle V \right \rangle ^{2}}{V \, k_{B} \, T}
222			\end{equation}
223
224			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.
225
226			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
227
228			\begin{equation}
229			\kappa_{T} = \frac{1}{N} \left ( \frac{\partial N}{\partial P} \right )_{T}
230			\end{equation}
231
232			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.
233
234			\subsection{Molecular orientation distribution at cluster boundary}
235
236			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.
237
238			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.
239
240			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.
241
242			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.
243
244			The orientation of a water molecule is described by
245
246			\begin{equation}
247	gezelter	3640	\cos{\theta}=\frac{\vec{r}_i\cdot\vec{\mu}_i}{\|\vec{r}_i\|\|\vec{\mu}_i\|}
248			\end{equation}
249
250	kstocke1	3649	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}.
251
252	gezelter	3640	\begin{figure}
253	kstocke1	3649	\includegraphics[width=\linewidth]{g_r_theta}
254			\caption{Definition of coordinates}
255			\label{coords}
256			\end{figure}
257
258			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).
259
260			\begin{figure}
261	gezelter	3640	\includegraphics[width=\linewidth]{pAngle}
262			\caption{SPC/E water clusters: only minor dewetting at the boundary}
263			\label{pAngle}
264			\end{figure}
265
266	kstocke1	3649	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.
267	gezelter	3640
268	kstocke1	3649	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.
269
270			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.
271
272
273	gezelter	3640	\subsection{Heterogeneous nanoparticle / water mixtures}
274
275
276			\section{Appendix A: Hydrodynamic tensor for triangular facets}
277
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
293			\section{Appendix B: Computing Convex Hulls on Parallel Computers}
294
295			\section{Acknowledgments}
296			Support for this project was provided by the
297			National Science Foundation under grant CHE-0848243. Computational
298			time was provided by the Center for Research Computing (CRC) at the
299			University of Notre Dame.
300
301			\newpage
302
303			\bibliography{langevinHull}
304
305			\end{doublespace}
306			\end{document}