--- trunk/oopsePaper/oopsePaper.tex	2004/07/28 15:44:21	1425
+++ trunk/oopsePaper/oopsePaper.tex	2004/07/29 18:01:05	1434
@@ -3,7 +3,7 @@
 \usepackage{amssymb}
 \usepackage{endfloat}
 \usepackage{listings}
-\usepackage{palatino}
+\usepackage{berkeley}
 \usepackage{graphicx}
 \usepackage[ref]{overcite}
 \usepackage{setspace}
@@ -35,21 +35,50 @@ Need an abstract.
 \maketitle
 
 \begin{abstract}
-Need an abstract.
+{\sc oopse} is a new molecular dynamics simulation program which is
+capable of efficiently integrating equations of motion for atom types
+with orientational degrees of freedom (e.g. ``sticky'' atoms and point
+dipoles).  Transition metals can also be simulated using the embedded
+atom method ({\sc eam}) potential included in the code.  Parallel
+simulations are carried out using the force-based decomposition
+method.  Simulations are specified using a very simple C-based
+meta-data language.  A number of advanced integrators are included,
+and the base integrator for orientational dynamics provides
+substantial improvements over older quaternion-base schemes.  All
+source code is available under a very permissive (BSD-style) Open
+Source license.
 \end{abstract}
 
 \section{\label{sec:intro}Introduction}
 
-UNDERWAY
+There are a number of excellent molecular dynamics packages available
+to the chemical physics
+community.\cite{Brooks83,MacKerell98,pearlman:1995,Gromacs,Gromacs3,DL_POLY,Tinker,Paradyn}
+All of these packages are stable, polished programs which solve many
+problems of interest.  Most are now capable of performing molecular
+dynamics simulations on parallel computers.  Some have source code
+which is freely available to the entire scientific community.  Few,
+however, are capable of efficiently integrating the equations of
+motion for atom types with orientational degrees of freedom
+(e.g. point dipoles, and ``sticky'' atoms).  And only one of the
+programs referenced can handle transition metal force fields like the
+Embedded Atom Method ({\sc eam}).  The direction our research program
+has taken us now involves the use of atoms with orientational degrees
+of freedom and transition metals.  Since these simulation methods may
+be of some use to other researchers, we have decided to release our
+program to the scientific community with a permissive open source
+license.
 
-
-We have structured this paper to first discuss the underlying concepts
-in this simulation package (Sec. \ref{oopseSec:IOfiles}).  The
-empirical energy functions implemented are discussed in
-Sec.~\ref{oopseSec:empiricalEnergy}.  Sec.~\ref{oopseSec:mechanics}
-describes the various Molecular Dynamics algorithms {\sc oopse}
-implements in the integration of the Newtonian equations of motion.
-Program design considerations for parallel computing are presented in
+This paper communicates the algorithmic details of our program, which
+we have been calling the Open source Object-oriented Parallel
+Simulation Engine (i.e. {\sc oopse}).  We have structured this paper
+to first discuss the underlying concepts in this simulation package
+(Sec. \ref{oopseSec:IOfiles}).  The empirical energy functions
+implemented are discussed in Sec.~\ref{oopseSec:empiricalEnergy}.
+Sec.~\ref{oopseSec:mechanics} describes the various Molecular Dynamics
+algorithms {\sc oopse} implements in the integration of Hamilton's
+equations of motion.  Program design considerations for parallel
+computing are presented in
 Sec.~\ref{oopseSec:parallelization}. Concluding remarks are presented
 in Sec.~\ref{oopseSec:conclusion}.
 
@@ -337,7 +366,7 @@ outlined in Sec.~\ref{oopseSec:empericalEnergy}. The f
 Sec.~\ref{oopseSec:minimizer}. The {\tt forceField} statement is
 important for the selection of which forces will be used in the course
 of the simulation. {\sc oopse} supports several force fields, as
-outlined in Sec.~\ref{oopseSec:empericalEnergy}. The force fields are
+outlined in Sec.~\ref{oopseSec:empiricalEnergy}. The force fields are
 interchangeable between simulations, with the only requirement being
 that all atoms needed by the simulation are defined within the
 selected force field. 
@@ -1048,6 +1077,35 @@ models can be found in reference~\citen{fennell04}. 
 and the drawbacks and benefits of the different density corrected SSD
 models can be found in reference~\citen{fennell04}. 
 
+\subsection{\label{oopseSec:WATER}The {\sc water} Force Field}
+
+In addition to the {\sc duff} force field's solvent description, a
+separate {\sc water} force field has been included for simulating many
+of the common rigid-body water models. In addition to the simple or
+dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD water), the common
+charge-based models were included (SPC, SPC/E, TIP3P, TIP4P, and
+TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00}
+In order to handle these models, charge-charge interactions were
+included in the force-loop:
+\begin{equation}
+V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, 
+\end{equation}
+where $q$ represents the charge on particle $i$ or $j$, and $e$ is the
+charge of an electron in Coulombs. The charge-charge interaction
+support is rudimentary in the current version of {\sc oopse}. As with
+the other pair interactions, charges can be simulated with a pure
+cutoff or a reaction field. The various methods for performing the
+Ewald summation have not yet been included. Also, the charge-dipole
+and charge-quadrupole (for interactions between SSD type water and
+charges) are not yet available, so it is currently inadvisable to mix
+dipolar and charge based molecules in the same system.
+
+The {\sc water} force field can be easily expanded through
+modification of the {\sc water} force field file ({\tt WATER.frc}). By
+adding atom types and inserting the appropriate parameters, it is
+possible to extend the force field to handle rigid molecules other
+than water.
+
 \subsection{\label{oopseSec:eam}Embedded Atom Method}
 
 {\sc oopse} implements a potential that describes bonding in
@@ -1409,37 +1467,75 @@ propagation. With the same time step, a 1000-molecule 
 
 The matrix rotations used in the DLM method end up being more costly
 computationally than the simpler arithmetic quaternion
-propagation. With the same time step, a 1000-molecule water simulation
-shows an average 7\% increase in computation time using the DLM method
-in place of quaternions. This cost is more than justified when
-comparing the energy conservation of the two methods as illustrated in
-Fig.~\ref{timestep}.
+propagation. With the same time step, a 1024-molecule water simulation
+shows an 12\% increase in computation time (averaged over several
+different time steps) using the DLM method in place of
+quaternions. This cost is more than justified when comparing the
+energy conservation of the two methods. Figure ~\ref{quatdlm} provides
+a comparative analysis of the {\sc dlm} method versus the simple
+quaternion method that was originally implemented.
 
 \begin{figure}
 \centering
-\includegraphics[width=\linewidth]{timeStep.eps}
-\caption[Energy conservation for quaternion versus DLM dynamics]{Energy conservation using quaternion based integration versus 
-the method proposed by Dullweber \emph{et al.} with increasing time
-step. For each time step, the dotted line is total energy using the
-DLM integrator, and the solid line comes from the quaternion
-integrator. The larger time step plots are shifted up from the true
-energy baseline for clarity.}
-\label{timestep}
+\includegraphics[width=\linewidth]{quatvsdlm.eps}
+\caption[Energy conservation analysis of the {\sc dlm} and quaternion 
+integration methods]{The logarithm of absolute value of the slope of
+the energy drift (\delta E$_1$) and the standard deviation of the
+energy fluctuations (\delta E$_0$) as a function of chosen time
+step. All simulations were of a 1024-molecule simulation of SSD water
+at 298 K starting from the same initial configuration. Note that the
+{\sc dlm} method provides a greater-than order-of-magnitude
+improvement in energy conservation and relative energy fluctuations
+over the quaternion method at all the tested time steps. The energy
+drift is quite steep for the larger time steps in both methods, and
+results in discontinuous behavior as the systems compound their
+anomolous energy accumulation.}
+\label{quatdlm}
 \end{figure}
 
-In Fig.~\ref{timestep}, the resulting energy drift at various time
-steps for both the DLM and quaternion integration schemes is
-compared. All of the 1000 molecule water simulations started with the
-same configuration, and the only difference was the method for
-handling rotational motion. At time steps of 0.1 and 0.5 fs, both
-methods for propagating molecule rotation conserve energy fairly well,
-with the quaternion method showing a slight energy drift over time in
-the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
-energy conservation benefits of the DLM method are clearly
-demonstrated. Thus, while maintaining the same degree of energy
-conservation, one can take considerably longer time steps, leading to
-an overall reduction in computation time.
+In Fig.~\ref{quatdlm}, \delta E$_1$ is a measure of the linear energy
+drift in units of kcal/mol per particle over a nanosecond of
+simulation time, and \delta E$_0$ is the standard deviation of the
+energy fluctuations in units of kcal/mol per particle. In the top
+plot, it is apparent that the energy drift is reduced by a significant
+amount (2 to 3 orders-of-magnitude improvement at every tested time
+step) by chosing the {\sc dlm} method over the simple non-symplectic
+quaternion integration method. When the energy drift becomes very
+small ($log_{10}[|\delta\text{E}_1|] < -3$), it is more difficult to
+calculate a slope, resulting in the larger displayed error bars. In
+addition to this improvement in energy drift, the fluctuation is the
+total energy are also dampened out by 1 to 2 orders-of-magnitude by
+utilizing the {\sc dlm} integration method.
 
+It was stated previously that the {\sc dlm} method was the more
+computationally expensive of the two implimented integration
+methodologies. In order to incorporate this information into the
+energy analysis a plot of energy drift versus computational cost was
+generated (Fig.~\ref{cpuCost}). This figure provides an estimate of
+the CPU time required under the two integration schemes for 1
+nanosecond of simulation time for the model 1024-molecule system. The
+plot is read by chosing a desired energy drift value and determining
+where both the curves cross. If a \delta E$_1$ of 1E-3 kcal/mol per
+particle is desired, a nanosecond of simulation time will require ~19
+hours of CPU time with the {\sc dlm} integrator, while the same small
+drift value will require ~154 hours of CPU time. This demonstrates the
+computational advantage of the integration scheme utilized in {\sc
+oopse}.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{compCost.eps}
+\caption[Energy drift as a function of required simulation run 
+time]{The logarithm of absolute value of the slope of the energy drift
+(\delta E$_1$) as a function of simulation run time. Simulations were
+performed on a single 2.5 GHz Pentium IV processor. Simulation time
+comparisons can be made by tracing horizontally from one curve to the
+other. For example, a simulation that takes ~24 hours using the {\sc
+dlm} method will take roughly 210 hours using the simple quaternion
+method if the same degree of energy conservation is desired.}
+\label{cpuCost}
+\end{figure}
+
 There is only one specific keyword relevant to the default integrator,
 and that is the time step for integrating the equations of motion.
 
@@ -2081,7 +2177,7 @@ used by default\\
 \end{table}
 
 
-\section{\label{sec:minimize}Energy Minimization}
+\section{\label{oopseSec:minimizer}Energy Minimization}
 
 As one of the basic procedures of molecular modeling, energy
 minimization is used to identify local configurations that are stable
