--- trunk/langevinHull/langevinHull.tex	2010/10/18 16:54:02	3651
+++ trunk/langevinHull/langevinHull.tex	2010/10/18 18:27:24	3652
@@ -41,17 +41,17 @@ Notre Dame, Indiana 46556}
   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.
+  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
@@ -76,10 +76,29 @@ Nos\'e-Hoover-Andersen equations of motion, the Berend
 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, the Berendsen pressure
-bath, and the Langevin Piston, all utilize coordinate transformation
-to adjust the box volume.
+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
+liquid which has a relatively uniform compressibility, the standard
+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 these approaches 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 collisions in the incompressible region.
+
 \begin{figure}
 \includegraphics[width=\linewidth]{AffineScale2}
 \caption{Affine Scaling constant pressure methods use box-length
@@ -92,9 +111,9 @@ to adjust the box volume.
 \label{affineScale}
 \end{figure}
 
+Additionally, one may often wish to simulate explicitly non-periodic
+systems, and the constraint that a periodic box must be used to 
 
-Heterogeneous mixtures of materials with different compressibilities?
-
 Explicitly non-periodic systems
 
 Elastic Bag 
@@ -103,71 +122,111 @@ A new method which uses a constant pressure and temper
 
 \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.
+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.
 
-Novel features: No a priori geometry is defined, No affine transforms,
-No fictitious particles, No bounding potentials.
+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.  
 
-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. 
+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
+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}
+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 a ($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}
-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)
+\Xi_f(t)\delta(t-t^\prime).
+\label{eq:randomForce}
 \end{eqnarray}
 
-Implemented in OpenMD.\cite{Meineke2005,openmd}
+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}).
 
+We have implemented this method by extending the Langevin dynamics
+integrator in our group code, OpenMD.\cite{Meineke2005,openmd}  
+
 \section{Tests \& Applications}
 
 \subsection{Bulk modulus of gold nanoparticles}