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