--- trunk/langevinHull/langevinHull.tex	2010/10/18 18:27:24	3652
+++ trunk/langevinHull/langevinHull.tex	2010/10/19 16:13:01	3653
@@ -80,15 +80,14 @@ transformation to adjust the box volume.
 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
+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.  
+atomic sites.
 
-The problem with these approaches becomes apparent when the material
+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
@@ -97,7 +96,7 @@ enough to avoid collisions in the incompressible regio
 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.
+enough to avoid the instabilities in the incompressible region.
 
 \begin{figure}
 \includegraphics[width=\linewidth]{AffineScale2}
@@ -111,16 +110,90 @@ Additionally, one may often wish to simulate explicitl
 \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 
+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]
 
-Explicitly non-periodic systems
+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]
 
-Elastic Bag 
+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?
 
-Spherical Boundary approaches
+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
@@ -191,11 +264,11 @@ and $\Xi_f(t)$ is a ($3 \times 3$) hydrodynamic tensor
 \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,
+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\
@@ -223,11 +296,56 @@ have the correct properties required by Eq. (\ref{eq:r
 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}  
+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}
 
@@ -334,21 +452,6 @@ The orientational preference exhibited by hull molecul
 
 \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}