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 solvent.
  We apply this method to several different systems including bare
  nanoparticles, nanoparticles in an 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{Melchionna1993} to the
Nos\'e-Hoover-Andersen equations of
motion,\cite{Hoover85,ANDERSEN:1980vn,Sturgeon:2000kx} the Berendsen
pressure bath,\cite{ISI:A1984TQ73500045} and the Langevin
Piston,\cite{FELLER:1995fk,Jakobsen:2005uq} all utilize coordinate
transformation to adjust the box volume.  As long as the material in
the simulation box is essentially a bulk-like liquid which has a
relatively uniform compressibility, the standard affine transform
approach provides an excellent way of adjusting the volume of the
system and applying pressure directly via the interactions between
atomic sites.

The problem with this approach becomes apparent when the material
being simulated is an inhomogeneous mixture in which portions of the
simulation box are incompressible relative to other portions.
Examples include simulations of metallic nanoparticles in liquid
environments, proteins at interfaces, as well as other multi-phase or
interfacial environments.  In these cases, the affine transform of
atomic coordinates will either cause numerical instability when the
sites in the incompressible medium collide with each other, or lead to
inefficient sampling of system volumes if the barostat is set slow
enough to avoid the instabilities in the incompressible region.

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

One may also wish to avoid affine transform periodic boundary methods
to simulate {\it explicitly non-periodic systems} under constant
pressure conditions. The use of periodic boxes to enforce a system
volume either requires effective solute concentrations that are much
higher than desirable, or unreasonable system sizes to avoid this
effect.  For example, calculations using typical hydration shells 
solvating a protein under periodic boundary conditions are quite
expensive. [CALCULATE EFFECTIVE PROTEIN CONCENTRATIONS IN TYPICAL
SIMULATIONS]

There have been a number of other approaches to explicit
non-periodicity that focus on constant or nearly-constant {\it volume}
conditions while maintaining bulk-like behavior.  Berkowitz and
McCammon introduced a stochastic (Langevin) boundary layer inside a
region of fixed molecules which effectively enforces constant
temperature and volume (NVT) conditions.\cite{Berkowitz1982} In this
approach, the stochastic and fixed regions were defined relative to a
central atom.  Brooks and Karplus extended this method to include
deformable stochastic boundaries.\cite{iii:6312} The stochastic
boundary approach has been used widely for protein
simulations. [CITATIONS NEEDED]

The electrostatic and dispersive behavior near the boundary has long
been a cause for concern.  King and Warshel introduced a surface
constrained all-atom solvent (SCAAS) which included polarization
effects of a fixed spherical boundary to mimic bulk-like behavior
without periodic boundaries.\cite{king:3647} In the SCAAS model, a
layer of fixed solvent molecules surrounds the solute and any explicit
solvent, and this in turn is surrounded by a continuum dielectric.
MORE HERE.  WHAT DID THEY FIND?

Beglov and Roux developed a boundary model in which the hard sphere
boundary has a radius that varies with the instantaneous configuration
of the solute (and solvent) molecules.\cite{beglov:9050} This model
contains a clear pressure and surface tension contribution to the free
energy which XXX.

Restraining {\it potentials} introduce repulsive potentials at the
surface of a sphere or other geometry.  The solute and any explicit
solvent are therefore restrained inside this potential.  Often the
potentials include a weak short-range attraction to maintain the
correct density at the boundary.  Beglov and Roux have also introduced
a restraining boundary potential which relaxes dynamically depending
on the solute geometry and the force the explicit system exerts on the
shell.\cite{Beglov:1995fk}

Recently, Krilov {\it et al.} introduced a flexible boundary model
that uses a Lennard-Jones potential between the solvent molecules and
a boundary which is determined dynamically from the position of the
nearest solute atom.\cite{LiY._jp046852t,Zhu:xw} This approach allows
the confining potential to prevent solvent molecules from migrating
too far from the solute surface, while providing a weak attractive
force pulling the solvent molecules towards a fictitious bulk solvent.
Although this approach is appealing and has physical motivation,
nanoparticles do not deform far from their original geometries even at
temperatures which vaporize the nearby solvent. For the systems like
the one described, the flexible boundary model will be nearly
identical to a fixed-volume restraining potential.

The approach of Kohanoff, Caro, and Finnis is the most promising of
the methods for introducing both constant pressure and temperature
into non-periodic simulations.\cite{Kohanoff:2005qm,Baltazar:2006ru}
This method is based on standard Langevin dynamics, but the Brownian
or random forces are allowed to act only on peripheral atoms and exert
force in a direction that is inward-facing relative to the facets of a
closed bounding surface.  The statistical distribution of the random
forces are uniquely tied to the pressure in the external reservoir, so
the method can be shown to sample the isobaric-isothermal ensemble.
Kohanoff {\it et al.} used a Delaunay tessellation to generate a
bounding surface surrounding the outermost atoms in the simulated
system.  This is not the only possible triangulated outer surface, but
guarantees that all of the random forces point inward towards the
cluster.

In the following sections, we extend and generalize the approach of
Kohanoff, Caro, and Finnis. The new method, which we are calling the
``Langevin Hull'' applies the external pressure, Langevin drag, and
random forces on the facets of the {\it hull itself} instead of the
atomic sites comprising the vertices of the hull.  This allows us to
decouple the external pressure contribution from the drag and random
force.  Section \ref{sec:meth}

\section{Methodology}
\label{sec:meth}

We have developed a new method which uses a constant pressure and
temperature bath.  This bath interacts only with the objects that are
currently at the edge of the system.  Since the edge is determined
dynamically as the simulation progresses, no {\it a priori} geometry
is defined.  The pressure and temperature bath interacts {\it
  directly} with the atoms on the edge and not with atoms interior to
the simulation.  This means that there are no affine transforms
required.  There are also no fictitious particles or bounding
potentials used in this approach.

The basics of the method are as follows. The simulation starts as a
collection of atomic locations in three dimensions (a point cloud).
Delaunay triangulation is used to find all facets between coplanar
neighbors.  In highly symmetric point clouds, facets can contain many
atoms, but in all but the most symmetric of cases one might experience
in a molecular dynamics simulation, the facets are simple triangles in
3-space that contain exactly three atoms.  

The convex hull is the set of facets that have {\it no concave
  corners} at an atomic site.  This eliminates all facets on the
interior of the point cloud, leaving only those exposed to the
bath. Sites on the convex hull are dynamic. As molecules re-enter the
cluster, all interactions between atoms on that molecule and the
external bath are removed.

For atomic sites in the interior of the point cloud, the equations of
motion are simple Newtonian dynamics,
\begin{equation}
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U,
\label{eq:Newton}
\end{equation}
where $m_i$ is the mass of site $i$, ${\mathbf v}_i(t)$ is the
instantaneous velocity of site $i$ at time $t$, and $U$ is the total
potential energy.  For atoms on the exterior of the cluster
(i.e. those that occupy one of the vertices of the convex hull), the
equation of motion is modified with an external force, ${\mathbf
  F}_i^{\mathrm ext}$,
\begin{equation}
m_i \dot{\mathbf v}_i(t)=-{\mathbf \nabla}_i U + {\mathbf F}_i^{\mathrm ext}.
\end{equation}

The external bath interacts directly with the facets of the convex
hull.  Since each vertex (or atom) provides one corner of a triangular
facet, the force on the facets are divided equally to each vertex.
However, each vertex can participate in multiple facets, so the resultant
force is a sum over all facets $f$ containing vertex $i$:
\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}

The external pressure bath applies a force to the facets of the convex
hull in direct proportion to the area of the facet, while the thermal
coupling depends on the solvent temperature, friction and the size and
shape of each facet. The thermal interactions are expressed as a
typical Langevin description of the forces,
\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}
Here, $P$ is the external pressure, $A_f$ and $\hat{n}_f$ are the area
and normal vectors for facet $f$, respectively.  ${\mathbf v}_f(t)$ is
the velocity of the facet,
\begin{equation}
{\mathbf v}_f(t) =  \frac{1}{3} \sum_{i=1}^{3} {\mathbf v}_i,
\end{equation}
and $\Xi_f(t)$ is an approximate ($3 \times 3$) hydrodynamic tensor
that depends on the geometry and surface area of facet $f$ and the
viscosity of the fluid (See Appendix A).  The hydrodynamic tensor is
related to the fluctuations of the random force, $\mathbf{R}(t)$, by
the fluctuation-dissipation theorem,
\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).
\label{eq:randomForce}
\end{eqnarray}

Once the hydrodynamic tensor is known for a given facet (see Appendix
A) obtaining a stochastic vector that has the properties in
Eq. (\ref{eq:randomForce}) can be done efficiently by carrying out a
one-time Cholesky decomposition to obtain the square root matrix of
the resistance tensor,
\begin{equation} 
\Xi_f = {\bf S} {\bf S}^{T},
\label{eq:Cholesky}
\end{equation} 
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
vector with the statistics required for the random force can then be
obtained by multiplying ${\bf S}$ onto a random 3-vector ${\bf Z}$ which
has elements chosen from a Gaussian distribution, such that:
\begin{equation}
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
\end{equation}
where $\delta t$ is the timestep in use during the simulation. The
random force, ${\bf R}_{f} = {\bf S} {\bf Z}$, can be shown to
have the correct properties required by Eq. (\ref{eq:randomForce}).

Our treatment of the hydrodynamic tensor must be approximate.  $\Xi$
for a triangular plate would normally be treated as a $6 \times 6$
tensor that includes translational and rotational drag as well as
translational-rotational coupling. The computation of hydrodynamic
tensors for rigid bodies has been detailed
elsewhere,\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002,Sun2008}
but the standard approach involving bead approximations would be
prohibitively expensive if it were recomputed at each step in a
molecular dynamics simulation.

We are utilizing an approximate hydrodynamic tensor obtained by first
constructing the Oseen tensor for the interaction of the centroid of
the facet ($f$) with each of the subfacets $j$,
\begin{equation}
T_{jf}=\frac{A_j}{8\pi\eta R_{jf}}\left(I +
  \frac{\mathbf{R}_{jf}\mathbf{R}_{jf}^T}{R_{jf}^2}\right)
\end{equation}
Here, $A_j$ is the area of subfacet $j$ which is a triangle containing
two of the vertices of the facet along with the centroid.
$\mathbf{R}_{jf}$ is the vector between the centroid of facet $f$ and
the centroid of sub-facet $j$, and $I$ is the ($3 \times 3$) identity
matrix.  $\eta$ is the viscosity of the external bath.

\begin{figure}
\includegraphics[width=\linewidth]{hydro}
\caption{The hydrodynamic tensor $\Xi$ for a facet comprising sites $i$,
  $j$, and $k$ is constructed using Oseen tensor contributions
  between the centoid of the facet $f$ and each of the sub-facets
  ($i,f,j$), ($j,f,k$), and ($k,f,i$). The centroids of the sub-facets
  are located at $1$, $2$, and $3$, and the area of each sub-facet is
  easily computed using half the cross product of two of the edges.}
\label{hydro}
\end{figure}

The Oseen tensors for each of the sub-facets are summed, and the
resulting matrix is inverted to give a $3 \times 3$ hydrodynamic
tensor for translations of the triangular plate,
\begin{equation}
\Xi_f(t) =\left[\sum_{i=1}^3 T_{if}\right]^{-1}.
\end{equation}
We have implemented this method by extending the Langevin dynamics
integrator in our group code, OpenMD.\cite{Meineke2005,openmd} There
is a moderate penalty for computing the convex hull at each step in
the molecular dynamics simulation (HOW MUCH?), but the convex hull is
remarkably easy to parallelize on distributed memory machines (see
Appendix B).

\section{Tests \& Applications}
\label{sec:tests}

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

The volume of a three-dimensional point cloud is not an obvious property to calculate. In order to calculate the isothermal compressibility we adapted the classic compressibility formula so that the compressibility could be calculated using information about the local density instead of the total volume of the convex hull. 

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


Assuming a uniform 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.

We initially used the classic compressibility formula 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{Debenedetti1986}, 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.


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

\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:	3655
Committed:	Tue Oct 19 18:25:01 2010 UTC (13 years, 8 months ago) by kstocke1
Content type:	application/x-tex
File size:	26599 byte(s)
Log Message:	M langevinHull.tex