| 22 |
|
simulation package of this size and scope would not have been possible |
| 23 |
|
without the collaborative efforts of my colleagues: Charles |
| 24 |
|
F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J.~Daniel |
| 25 |
< |
Gezelter. Although my contributions to [\sc oopse} are signifigant, |
| 25 |
> |
Gezelter. Although my contributions to {\sc oopse} are significant, |
| 26 |
|
consideration of my work apart from the others, would not give a |
| 27 |
|
complete description to the package's capabilities. As such, all |
| 28 |
|
contributions to {\sc oopse} to date are presented in this chapter. |
| 29 |
|
|
| 30 |
< |
{\sc give final breakdown of who wrote which section here.} |
| 30 |
> |
Charles Vardeman is responsible for the parallelization of {\sc oopse} |
| 31 |
> |
(Sec.~\ref{oopseSec:parallelization}) as well as the inclusion of the |
| 32 |
> |
embedded-atom potential (Sec.~\ref{oopseSec:eam}). Teng Lin's |
| 33 |
> |
contributions include refinement of the periodic boundary conditions |
| 34 |
> |
(Sec.~\ref{oopseSec:pbc}), the z-constraint method |
| 35 |
> |
(Sec.~\ref{oopseSec:zcons}), refinement of the property analysis |
| 36 |
> |
programs (Sec.~\ref{oopseSec:props}), and development in the extended |
| 37 |
> |
state integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher |
| 38 |
> |
Fennell worked on the symplectic integrator |
| 39 |
> |
(Sec.~\ref{oopseSec:integrate}) and the refinement of the {\sc ssd} |
| 40 |
> |
water model (Sec.~\ref{oopseSec:SSD}). Daniel Gezelter lent his |
| 41 |
> |
talents in the development of the extended system integrators |
| 42 |
> |
(Sec.~\ref{oopseSec:noseHooverThermo}) as well as giving general |
| 43 |
> |
direction and oversight to the entire project. My responsibilities |
| 44 |
> |
covered the creation and specification of {\sc bass} |
| 45 |
> |
(Sec.~\ref{oopseSec:IOfiles}), the original development of the single |
| 46 |
> |
processor version of {\sc oopse}, contributions to the extended state |
| 47 |
> |
integrators (Sec.~\ref{oopseSec:noseHooverThermo}), the implementation |
| 48 |
> |
of the Lennard-Jones (Sec.~\ref{sec:LJPot}) and {\sc duff} |
| 49 |
> |
(Sec.~\ref{oopseSec:DUFF}) force fields, and initial implementation of |
| 50 |
> |
the property analysis (Sec.~\ref{oopseSec:props}) and system |
| 51 |
> |
initialization (Sec.~\ref{oopseSec:initCoords}) utility programs. |
| 52 |
|
|
| 53 |
|
\section{\label{sec:intro}Introduction} |
| 54 |
|
|
| 57 |
|
might consider writing one's own programming code. However, as systems |
| 58 |
|
grow larger and more complex, building and maintaining code for the |
| 59 |
|
simulations becomes a time consuming task. In such cases it is usually |
| 60 |
< |
more convienent for a researcher to turn to pre-existing simulation |
| 60 |
> |
more convenient for a researcher to turn to pre-existing simulation |
| 61 |
|
packages. These packages, such as {\sc amber}\cite{pearlman:1995} and |
| 62 |
|
{\sc charmm}\cite{Brooks83}, provide powerful tools for researchers to |
| 63 |
|
conduct simulations of their systems without spending their time |
| 64 |
|
developing a code base to conduct their research. This then frees them |
| 65 |
< |
to perhaps explore experimental analouges to their models. |
| 65 |
> |
to perhaps explore experimental analogues to their models. |
| 66 |
|
|
| 67 |
|
Despite their utility, problems with these packages arise when |
| 68 |
|
researchers try to develop techniques or energetic models that the |
| 69 |
|
code was not originally designed to do. Examples of uncommonly |
| 70 |
|
implemented techniques and energetics include; dipole-dipole |
| 71 |
< |
interactions, rigid body dynamics, and metallic emmbedded |
| 71 |
> |
interactions, rigid body dynamics, and metallic embedded |
| 72 |
|
potentials. When faced with these obstacles, a researcher must either |
| 73 |
|
develop their own code or license and extend one of the commercial |
| 74 |
|
packages. What we have elected to do, is develop a package of |
| 76 |
|
our research is based. |
| 77 |
|
|
| 78 |
|
Having written {\sc oopse} we are implementing the concept of Open |
| 79 |
< |
Source dcevelopment, and releaseing our source code into the public |
| 79 |
> |
Source development, and releasing our source code into the public |
| 80 |
|
domain. It is our intent that by doing so, other researchers might |
| 81 |
|
benefit from our work, and add their own contributions to the |
| 82 |
|
package. The license under which {\sc oopse} is distributed allows any |
| 85 |
|
only the models of interest to ourselves, but also those of the |
| 86 |
|
community of scientists who contribute back to the project. |
| 87 |
|
|
| 88 |
< |
We have structured this chapter to first discuss the emperical energy |
| 88 |
> |
We have structured this chapter to first discuss the empirical energy |
| 89 |
|
functions that {\sc oopse } implements in |
| 90 |
< |
Sec.~\ref{oopseSec:empericalEnergy}. Following that is a discusion of |
| 90 |
> |
Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of |
| 91 |
|
the various input and output files associated with the package |
| 92 |
< |
(Sec.~\ref{oopseSec:IOfiles}). In Sec.~\ref{oopseSec:Mechanics} |
| 92 |
> |
(Sec.~\ref{oopseSec:IOfiles}). In Sec.~\ref{oopseSec:mechanics} |
| 93 |
|
elucidates the various Molecular Dynamics algorithms {\sc oopse} |
| 94 |
< |
mplements in the integration of the Newtonian equations of |
| 94 |
> |
implements in the integration of the Newtonian equations of |
| 95 |
|
motion. Basic analysis of the trajectories obtained from the |
| 96 |
|
simulation is discussed in Sec.~\ref{oopseSec:props}. Program design |
| 97 |
|
considerations as well as the software distribution license is |
| 108 |
|
atoms of an element, or be used for collections of atoms such as |
| 109 |
|
methyl and carbonyl groups. The atoms are also capable of having |
| 110 |
|
directional components associated with them (\emph{e.g.}~permanent |
| 111 |
< |
dipoles). Charges on atoms are not currently supported by {\sc oopse}. |
| 111 |
> |
dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for |
| 112 |
> |
a given atom type are set in the force field parameter files |
| 113 |
|
|
| 114 |
< |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of the simple Ar molecule},label=sch:AtmMole] |
| 114 |
> |
\begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole] |
| 115 |
|
molecule{ |
| 116 |
|
name = "Ar"; |
| 117 |
|
nAtoms = 1; |
| 123 |
|
\end{lstlisting} |
| 124 |
|
|
| 125 |
|
|
| 126 |
< |
Atoms can be collected into secondary srtructures such as rigid bodies |
| 126 |
> |
Atoms can be collected into secondary structures such as rigid bodies |
| 127 |
|
or molecules. The molecule is a way for {\sc oopse} to keep track of |
| 128 |
|
the atoms in a simulation in logical manner. Molecular units store the |
| 129 |
< |
identities of all the atoms associated with themselves, and are |
| 130 |
< |
responsible for the evaluation of their own internal interactions |
| 131 |
< |
(\emph{i.e.}~bonds, bends, and torsions). Scheme \ref{sch:AtmMole} |
| 132 |
< |
shws how one creates a molecule in the \texttt{.mdl} files. The |
| 133 |
< |
position of the atoms given in the declaration are relative to the |
| 134 |
< |
origin of the molecule, and is used when creating a system containing |
| 135 |
< |
the molecule. |
| 129 |
> |
identities of all the atoms and rigid bodies associated with |
| 130 |
> |
themselves, and are responsible for the evaluation of their own |
| 131 |
> |
internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme |
| 132 |
> |
\ref{sch:AtmMole} shows how one creates a molecule in a |
| 133 |
> |
\texttt{.mdl} file. The position of the atoms given in the |
| 134 |
> |
declaration are relative to the origin of the molecule, and is used |
| 135 |
> |
when creating a system containing the molecule. |
| 136 |
|
|
| 137 |
|
As stated previously, one of the features that sets {\sc oopse} apart |
| 138 |
|
from most of the current molecular simulation packages is the ability |
| 180 |
|
|
| 181 |
|
{\sc oopse} utilizes a relatively new scheme that propagates the |
| 182 |
|
entire nine parameter rotation matrix internally. Further discussion |
| 183 |
< |
on this choice can be found in Sec.~\ref{sec:integrate}. An example |
| 184 |
< |
definition of a riged body can be seen in Scheme |
| 183 |
> |
on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example |
| 184 |
> |
definition of a rigid body can be seen in Scheme |
| 185 |
|
\ref{sch:rigidBody}. The positions in the atom definitions are the |
| 186 |
|
placements of the atoms relative to the origin of the rigid body, |
| 187 |
|
which itself has a position relative to the origin of the molecule. |
| 210 |
|
} |
| 211 |
|
\end{lstlisting} |
| 212 |
|
|
| 213 |
< |
\subsection{\label{sec:LJPot}The Lennard Jones Potential} |
| 213 |
> |
\subsection{\label{sec:LJPot}The Lennard Jones Force Field} |
| 214 |
|
|
| 215 |
|
The most basic force field implemented in {\sc oopse} is the |
| 216 |
< |
Lennard-Jones potential, which mimics the van der Waals interaction at |
| 216 |
> |
Lennard-Jones force field, which mimics the van der Waals interaction at |
| 217 |
|
long distances, and uses an empirical repulsion at short |
| 218 |
|
distances. The Lennard-Jones potential is given by: |
| 219 |
|
\begin{equation} |
| 227 |
|
Where $r_{ij}$ is the distance between particles $i$ and $j$, |
| 228 |
|
$\sigma_{ij}$ scales the length of the interaction, and |
| 229 |
|
$\epsilon_{ij}$ scales the well depth of the potential. Scheme |
| 230 |
< |
\ref{sch:LJFF} gives and example partial \texttt{.bass} file that |
| 231 |
< |
shows a system of 108 Ar particles simulated with the Lennard-Jones |
| 232 |
< |
force field. |
| 230 |
> |
\ref{sch:LJFF} gives and example \texttt{.bass} file that |
| 231 |
> |
sets up a system of 108 Ar particles to be simulated using the |
| 232 |
> |
Lennard-Jones force field. |
| 233 |
|
|
| 234 |
|
\begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}] |
| 235 |
|
|
| 288 |
|
The dipolar unified-atom force field ({\sc duff}) was developed to |
| 289 |
|
simulate lipid bilayers. The simulations require a model capable of |
| 290 |
|
forming bilayers, while still being sufficiently computationally |
| 291 |
< |
efficient to allow large systems ($\approx$100's of phospholipids, |
| 292 |
< |
$\approx$1000's of waters) to be simulated for long times |
| 293 |
< |
($\approx$10's of nanoseconds). |
| 291 |
> |
efficient to allow large systems ($\sim$100's of phospholipids, |
| 292 |
> |
$\sim$1000's of waters) to be simulated for long times |
| 293 |
> |
($\sim$10's of nanoseconds). |
| 294 |
|
|
| 295 |
|
With this goal in mind, {\sc duff} has no point |
| 296 |
|
charges. Charge-neutral distributions were replaced with dipoles, |
| 306 |
|
the head group center of mass, our model mimics the head group of |
| 307 |
|
phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones |
| 308 |
|
site is located at the pseudoatom's center of mass. The model is |
| 309 |
< |
illustrated by the dark grey atom in Fig.~\ref{fig:lipidModel}. The |
| 309 |
> |
illustrated by the dark grey atom in Fig.~\ref{oopseFig:lipidModel}. The |
| 310 |
|
water model we use to complement the dipoles of the lipids is our |
| 311 |
|
reparameterization of the soft sticky dipole (SSD) model of Ichiye |
| 312 |
|
\emph{et al.}\cite{liu96:new_model} |
| 367 |
|
The total potential energy function in {\sc duff} is |
| 368 |
|
\begin{equation} |
| 369 |
|
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
| 370 |
< |
+ \sum^{N}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} |
| 370 |
> |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} |
| 371 |
|
\label{eq:totalPotential} |
| 372 |
|
\end{equation} |
| 373 |
|
Where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
| 394 |
|
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot} |
| 395 |
|
\end{equation} |
| 396 |
|
Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
| 397 |
< |
(see Fig.~\ref{fig:lipidModel}), $\theta_0$ is the equilibrium |
| 397 |
> |
(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium |
| 398 |
|
bond angle, and $k_{\theta}$ is the force constant which determines the |
| 399 |
|
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
| 400 |
|
$\theta_0$ are borrowed from those in TraPPE.\cite{Siepmann1998} |
| 408 |
|
\label{eq:origTorsionPot} |
| 409 |
|
\end{equation} |
| 410 |
|
Here $\phi$ is the angle defined by four bonded neighbors $i$, |
| 411 |
< |
$j$, $k$, and $l$ (again, see Fig.~\ref{fig:lipidModel}). For |
| 411 |
> |
$j$, $k$, and $l$ (again, see Fig.~\ref{oopseFig:lipidModel}). For |
| 412 |
|
computational efficiency, the torsion potential has been recast after |
| 413 |
< |
the method of CHARMM,\cite{charmm1983} in which the angle series is |
| 413 |
> |
the method of {\sc charmm},\cite{Brooks83} in which the angle series is |
| 414 |
|
converted to a power series of the form: |
| 415 |
|
\begin{equation} |
| 416 |
|
V_{\text{torsion}}(\phi) = |
| 443 |
|
Where $V_{\text{LJ}}$ is the Lennard Jones potential, |
| 444 |
|
$V_{\text{dipole}}$ is the dipole dipole potential, and |
| 445 |
|
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
| 446 |
< |
(Sec.~\ref{sec:SSD}). Note that not all atom types include all |
| 446 |
> |
(Sec.~\ref{oopseSec:SSD}). Note that not all atom types include all |
| 447 |
|
interactions. |
| 448 |
|
|
| 449 |
|
The dipole-dipole potential has the following form: |
| 461 |
|
towards $j$, and $\boldsymbol{\Omega}_i$ and $\boldsymbol{\Omega}_j$ |
| 462 |
|
are the orientational degrees of freedom for atoms $i$ and $j$ |
| 463 |
|
respectively. $|\mu_i|$ is the magnitude of the dipole moment of atom |
| 464 |
< |
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation |
| 465 |
< |
vector of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is |
| 466 |
< |
the unit vector pointing along $\mathbf{r}_{ij}$. |
| 464 |
> |
$i$, $\boldsymbol{\hat{u}}_i$ is the standard unit orientation vector |
| 465 |
> |
of $\boldsymbol{\Omega}_i$, and $\boldsymbol{\hat{r}}_{ij}$ is the |
| 466 |
> |
unit vector pointing along $\mathbf{r}_{ij}$ |
| 467 |
> |
($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$). |
| 468 |
|
|
| 469 |
|
|
| 470 |
< |
\subsubsection{\label{sec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 470 |
> |
\subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF} |
| 471 |
|
|
| 472 |
|
In the interest of computational efficiency, the default solvent used |
| 473 |
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
| 553 |
|
model. Solvent parameters can be easily modified in an accompanying |
| 554 |
|
{\sc BASS} file as illustrated in the scheme below. A table of the |
| 555 |
|
parameter values and the drawbacks and benefits of the different |
| 556 |
< |
density corrected SSD models can be found in reference |
| 557 |
< |
\ref{Gezelter04}. |
| 556 |
> |
density corrected SSD models can be found in |
| 557 |
> |
reference~\cite{Gezelter04}. |
| 558 |
|
|
| 559 |
|
\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]An example file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
| 560 |
|
|
| 593 |
|
|
| 594 |
|
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
| 595 |
|
|
| 596 |
< |
Several other molecular dynamics packages\cite{dynamo86} exist which have the |
| 596 |
> |
There are Molecular Dynamics packages which have the |
| 597 |
|
capacity to simulate metallic systems, including some that have |
| 598 |
|
parallel computational abilities\cite{plimpton93}. Potentials that |
| 599 |
|
describe bonding transition metal |
| 621 |
|
spherically averaged atomic electron densities given by |
| 622 |
|
$\rho_{i}$. $\phi_{ij}$ is a primarily repulsive pairwise interaction |
| 623 |
|
between atoms $i$ and $j$. In the original formulation of |
| 624 |
< |
{\sc eam} cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however |
| 624 |
> |
{\sc eam}\cite{Daw84}, $\phi_{ij}$ was an entirely repulsive term, however |
| 625 |
|
in later refinements to EAM have shown that non-uniqueness between $F$ |
| 626 |
|
and $\phi$ allow for more general forms for $\phi$.\cite{Daw89} |
| 627 |
|
There is a cutoff distance, $r_{cut}$, which limits the |
| 628 |
|
summations in the {\sc eam} equation to the few dozen atoms |
| 629 |
|
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
| 630 |
< |
interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FDB86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}. |
| 630 |
> |
interactions. Foiles et al. fit EAM potentials for fcc metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals\cite{FBD86}. These potential fits are in the DYNAMO 86 format and are included with {\sc oopse}. |
| 631 |
|
|
| 632 |
|
|
| 633 |
|
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
| 694 |
|
|
| 695 |
|
\subsection{{\sc bass} and Model Files} |
| 696 |
|
|
| 697 |
< |
Every {\sc oopse} simuation begins with a {\sc bass} file. {\sc bass} |
| 697 |
> |
Every {\sc oopse} simulation begins with a {\sc bass} file. {\sc bass} |
| 698 |
|
(\underline{B}izarre \underline{A}tom \underline{S}imulation |
| 699 |
|
\underline{S}yntax) is a script syntax that is parsed by {\sc oopse} at |
| 700 |
|
runtime. The {\sc bass} file allows for the user to completely describe the |
| 710 |
|
\label{fig:bassExample} |
| 711 |
|
\end{figure} |
| 712 |
|
|
| 713 |
< |
Within the \texttt{.bass} file it is neccassary to provide a complete |
| 713 |
> |
Within the \texttt{.bass} file it is necessary to provide a complete |
| 714 |
|
description of the molecule before it is actually placed in the |
| 715 |
|
simulation. The {\sc bass} syntax was originally developed with this goal in |
| 716 |
|
mind, and allows for the specification of all the atoms in a molecular |
| 717 |
|
prototype, as well as any bonds, bends, or torsions. These |
| 718 |
|
descriptions can become lengthy for complex molecules, and it would be |
| 719 |
< |
inconvient to duplicate the simulation at the begining of each {\sc bass} |
| 719 |
> |
inconvenient to duplicate the simulation at the beginning of each {\sc bass} |
| 720 |
|
script. Addressing this issue {\sc bass} allows for the inclusion of model |
| 721 |
|
files at the top of a \texttt{.bass} file. These model files, denoted |
| 722 |
|
with the \texttt{.mdl} extension, allow the user to describe a |
| 735 |
|
There are three major files used by {\sc oopse} written in the coordinate |
| 736 |
|
format, they are as follows: the initialization file, the simulation |
| 737 |
|
trajectory file, and the final coordinates of the simulation. The |
| 738 |
< |
initialization file is neccassary for {\sc oopse} to start the simulation |
| 738 |
> |
initialization file is necessary for {\sc oopse} to start the simulation |
| 739 |
|
with the proper coordinates. It is typically denoted with the |
| 740 |
|
extension \texttt{.init}. The trajectory file is created at the |
| 741 |
|
beginning of the simulation, and is used to store snapshots of the |
| 744 |
|
file at an interval specified in the \texttt{.bass} file. The |
| 745 |
|
trajectory file is given the extension \texttt{.dump}. The final |
| 746 |
|
coordinate file is the end of run or \texttt{.eor} file. The |
| 747 |
< |
\texttt{.eor} file stores the final configuration of teh system for a |
| 747 |
> |
\texttt{.eor} file stores the final configuration of the system for a |
| 748 |
|
given simulation. The file is updated at the same time as the |
| 749 |
|
\texttt{.dump} file. However, it only contains the most recent |
| 750 |
|
frame. In this way, an \texttt{.eor} file may be used as the |
| 751 |
|
initialization file to a second simulation in order to continue or |
| 752 |
|
recover the previous simulation. |
| 753 |
|
|
| 754 |
< |
\subsection{Generation of Initial Coordinates} |
| 754 |
> |
\subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates} |
| 755 |
|
|
| 756 |
< |
As was stated in Sec.~\ref{subSec:coordFiles}, an initialization file |
| 756 |
> |
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization file |
| 757 |
|
is needed to provide the starting coordinates for a simulation. The |
| 758 |
|
{\sc oopse} package provides a program called \texttt{sysBuilder} to aid in |
| 759 |
|
the creation of the \texttt{.init} file. \texttt{sysBuilder} is {\sc bass} |
| 760 |
|
aware, and will recognize arguments and parameters in the |
| 761 |
|
\texttt{.bass} file that would otherwise be ignored by the |
| 762 |
< |
simulation. The program itself is under contiunual development, and is |
| 762 |
> |
simulation. The program itself is under continual development, and is |
| 763 |
|
offered here as a helper tool only. |
| 764 |
|
|
| 765 |
|
\subsection{The Statistics File} |
| 768 |
|
file records such statistical quantities as the instantaneous |
| 769 |
|
temperature, volume, pressure, etc. It is written out with the |
| 770 |
|
frequency specified in the \texttt{.bass} file. The file allows the |
| 771 |
< |
user to observe the system variables as a function od simulation time |
| 771 |
> |
user to observe the system variables as a function of simulation time |
| 772 |
|
while the simulation is in progress. One useful function the |
| 773 |
|
statistics file serves is to monitor the conserved quantity of a given |
| 774 |
|
simulation ensemble, this allows the user to observe the stability of |
| 777 |
|
|
| 778 |
|
\section{\label{oopseSec:mechanics}Mechanics} |
| 779 |
|
|
| 780 |
< |
\subsection{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 780 |
> |
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 781 |
|
|
| 782 |
|
Integration of the equations of motion was carried out using the |
| 783 |
|
symplectic splitting method proposed by Dullweber \emph{et |
| 856 |
|
{\sc oopse} implements a |
| 857 |
|
|
| 858 |
|
|
| 859 |
< |
\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 859 |
> |
\subsubsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 860 |
|
|
| 861 |
|
To mimic the effects of being in a constant temperature ({\sc nvt}) |
| 862 |
|
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
| 884 |
|
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
| 885 |
|
|
| 886 |
|
|
| 887 |
< |
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 887 |
> |
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 888 |
|
|
| 889 |
< |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 889 |
> |
Based on fluctuation-dissipation theorem,\bigskip\ force auto-correlation |
| 890 |
|
method was developed to investigate the dynamics of ions inside the ion |
| 891 |
|
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 892 |
< |
from the deviation of the instaneous force from its mean force. |
| 892 |
> |
from the deviation of the instantaneous force from its mean force. |
| 893 |
|
|
| 894 |
|
% |
| 895 |
|
|
| 925 |
|
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 926 |
|
well as subtract the total constraint forces from the rest of the system after |
| 927 |
|
force calculation at each time step. |
| 928 |
< |
\begin{verbatim} |
| 929 |
< |
$F_{\alpha i}=0$ |
| 930 |
< |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$ |
| 931 |
< |
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$ |
| 932 |
< |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}$ |
| 933 |
< |
\end{verbatim} |
| 928 |
> |
\begin{align} |
| 929 |
> |
F_{\alpha i}&=0\\ |
| 930 |
> |
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}\\ |
| 931 |
> |
F_{\alpha i}&=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}\\ |
| 932 |
> |
V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}} |
| 933 |
> |
\end{align} |
| 934 |
|
|
| 935 |
|
At the very beginning of the simulation, the molecules may not be at its |
| 936 |
|
constraint position. To move the z-constraint molecule to the specified |
| 1080 |
|
between the frames contained within the two blocks are |
| 1081 |
|
calculated. Once completed, the memory blocks are incremented, and the |
| 1082 |
|
process is repeated. A diagram illustrating the process is shown in |
| 1083 |
< |
Fig.~\ref{fig:dynamicPropsMemory}. As was the case with \texttt{staticProps}, |
| 1084 |
< |
multiple properties may be calculated in a single run to avoid |
| 1085 |
< |
multiple reads on the same file. |
| 1083 |
> |
Fig.~\ref{oopseFig:dynamicPropsMemory}. As was the case with |
| 1084 |
> |
\texttt{staticProps}, multiple properties may be calculated in a |
| 1085 |
> |
single run to avoid multiple reads on the same file. |
| 1086 |
|
|
| 1064 |
– |
\begin{figure} |
| 1065 |
– |
\centering |
| 1066 |
– |
\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 1067 |
– |
\caption{This diagram illustrates the dynamic memory allocation used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.} |
| 1068 |
– |
\label{fig:dynamicPropsMemory} |
| 1069 |
– |
\end{figure} |
| 1087 |
|
|
| 1088 |
+ |
|
| 1089 |
|
\section{\label{oopseSec:design}Program Design} |
| 1090 |
|
|
| 1091 |
< |
\subsection{\label{sec:architecture} OOPSE Architecture} |
| 1091 |
> |
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
| 1092 |
|
|
| 1093 |
< |
The core of OOPSE is divided into two main object libraries: {\texttt |
| 1094 |
< |
libBASS} and {\texttt libmdtools}. {\texttt libBASS} is the library |
| 1095 |
< |
developed around the parseing engine and {\texttt libmdtools} is the |
| 1096 |
< |
software library developed around the simulation engine. |
| 1093 |
> |
The core of OOPSE is divided into two main object libraries: |
| 1094 |
> |
\texttt{libBASS} and \texttt{libmdtools}. \texttt{libBASS} is the |
| 1095 |
> |
library developed around the parsing engine and \texttt{libmdtools} |
| 1096 |
> |
is the software library developed around the simulation engine. These |
| 1097 |
> |
two libraries are designed to encompass all the basic functions and |
| 1098 |
> |
tools that {\sc oopse} provides. Utility programs, such as the |
| 1099 |
> |
property analyzers, need only link against the software libraries to |
| 1100 |
> |
gain access to parsing, force evaluation, and input / output |
| 1101 |
> |
routines. |
| 1102 |
|
|
| 1103 |
< |
\subsection{\label{sec:parallelization} Parallelization of OOPSE} |
| 1103 |
> |
Contained in \texttt{libBASS} are all the routines associated with |
| 1104 |
> |
reading and parsing the \texttt{.bass} input files. Given a |
| 1105 |
> |
\texttt{.bass} file, \texttt{libBASS} will open it and any associated |
| 1106 |
> |
\texttt{.mdl} files; then create structures in memory that are |
| 1107 |
> |
templates of all the molecules specified in the input files. In |
| 1108 |
> |
addition, any simulation parameters set in the \texttt{.bass} file |
| 1109 |
> |
will be placed in a structure for later query by the controlling |
| 1110 |
> |
program. |
| 1111 |
|
|
| 1112 |
< |
Although processor power is doubling roughly every 18 months according |
| 1113 |
< |
to the famous Moore's Law\cite{moore}, it is still unreasonable to |
| 1114 |
< |
simulate systems of more then a 1000 atoms on a single processor. To |
| 1115 |
< |
facilitate study of larger system sizes or smaller systems on long |
| 1116 |
< |
time scales in a reasonable period of time, parallel methods were |
| 1117 |
< |
developed allowing multiple CPU's to share the simulation |
| 1118 |
< |
workload. Three general categories of parallel decomposition method's |
| 1119 |
< |
have been developed including atomic, spatial and force decomposition |
| 1120 |
< |
methods. |
| 1112 |
> |
Located in \texttt{libmdtools} are all other routines necessary to a |
| 1113 |
> |
Molecular Dynamics simulation. The library uses the main data |
| 1114 |
> |
structures returned by \texttt{libBASS} to initialize the various |
| 1115 |
> |
parts of the simulation: the atom structures and positions, the force |
| 1116 |
> |
field, the integrator, \emph{et cetera}. After initialization, the |
| 1117 |
> |
library can be used to perform a variety of tasks: integrate a |
| 1118 |
> |
Molecular Dynamics trajectory, query phase space information from a |
| 1119 |
> |
specific frame of a completed trajectory, or even recalculate force or |
| 1120 |
> |
energetic information about specific frames from a completed |
| 1121 |
> |
trajectory. |
| 1122 |
|
|
| 1123 |
+ |
With these core libraries in place, several programs have been |
| 1124 |
+ |
developed to utilize the routines provided by \texttt{libBASS} and |
| 1125 |
+ |
\texttt{libmdtools}. The main program of the package is \texttt{oopse} |
| 1126 |
+ |
and the corresponding parallel version \texttt{oopse\_MPI}. These two |
| 1127 |
+ |
programs will take the \texttt{.bass} file, and create then integrate |
| 1128 |
+ |
the simulation specified in the script. The two analysis programs |
| 1129 |
+ |
\texttt{staticProps} and \texttt{dynamicProps} utilize the core |
| 1130 |
+ |
libraries to initialize and read in trajectories from previously |
| 1131 |
+ |
completed simulations, in addition to the ability to use functionality |
| 1132 |
+ |
from \texttt{libmdtools} to recalculate forces and energies at key |
| 1133 |
+ |
frames in the trajectories. Lastly, the family of system building |
| 1134 |
+ |
programs (Sec.~\ref{oopseSec:initCoords}) also use the libraries to |
| 1135 |
+ |
store and output the system configurations they create. |
| 1136 |
+ |
|
| 1137 |
+ |
\subsection{\label{oopseSec:parallelization} Parallelization of {\sc oopse}} |
| 1138 |
+ |
|
| 1139 |
+ |
Although processor power is continually growing month by month, it is |
| 1140 |
+ |
still unreasonable to simulate systems of more then a 1000 atoms on a |
| 1141 |
+ |
single processor. To facilitate study of larger system sizes or |
| 1142 |
+ |
smaller systems on long time scales in a reasonable period of time, |
| 1143 |
+ |
parallel methods were developed allowing multiple CPU's to share the |
| 1144 |
+ |
simulation workload. Three general categories of parallel |
| 1145 |
+ |
decomposition method's have been developed including atomic, spatial |
| 1146 |
+ |
and force decomposition methods. |
| 1147 |
+ |
|
| 1148 |
|
Algorithmically simplest of the three method's is atomic decomposition |
| 1149 |
|
where N particles in a simulation are split among P processors for the |
| 1150 |
|
duration of the simulation. Computational cost scales as an optimal |
| 1179 |
|
favorably then spatial decomposition up to 10,000 atoms and favorably |
| 1180 |
|
competes with spatial methods for up to 100,000 atoms. |
| 1181 |
|
|
| 1182 |
< |
\subsection{\label{openSource}Open Source and Distribution License} |
| 1182 |
> |
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 1183 |
|
|
| 1184 |
+ |
For large simulations, the trajectory files can sometimes reach sizes |
| 1185 |
+ |
in excess of several gigabytes. In order to effectively analyze that |
| 1186 |
+ |
amount of data+, two memory management schemes have been devised for |
| 1187 |
+ |
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
| 1188 |
+ |
developed for \texttt{staticProps}, is the simplest. As each frame's |
| 1189 |
+ |
statistics are calculated independent of each other, memory is |
| 1190 |
+ |
allocated for each frame, then freed once correlation calculations are |
| 1191 |
+ |
complete for the snapshot. To prevent multiple passes through a |
| 1192 |
+ |
potentially large file, \texttt{staticProps} is capable of calculating |
| 1193 |
+ |
all requested correlations per frame with only a single pair loop in |
| 1194 |
+ |
each frame and a single read through of the file. |
| 1195 |
|
|
| 1196 |
< |
\section{\label{oopseSec:conclusion}Conclusion} |
| 1196 |
> |
The second, more advanced memory scheme, is used by |
| 1197 |
> |
\texttt{dynamicProps}. Here, the program must have multiple frames in |
| 1198 |
> |
memory to calculate time dependent correlations. In order to prevent a |
| 1199 |
> |
situation where the program runs out of memory due to large |
| 1200 |
> |
trajectories, the user is able to specify that the trajectory be read |
| 1201 |
> |
in blocks. The number of frames in each block is specified by the |
| 1202 |
> |
user, and upon reading a block of the trajectory, |
| 1203 |
> |
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 1204 |
> |
pairs within the block. After in block correlations are complete, a |
| 1205 |
> |
second block of the trajectory is read, and the cross correlations are |
| 1206 |
> |
calculated between the two blocks. this second block is then freed and |
| 1207 |
> |
then incremented and the process repeated until the end of the |
| 1208 |
> |
trajectory. Once the end is reached, the first block is freed then |
| 1209 |
> |
incremented, and the again the internal time correlations are |
| 1210 |
> |
calculated. The algorithm with the second block is then repeated with |
| 1211 |
> |
the new origin block, until all frame pairs have been correlated in |
| 1212 |
> |
time. This process is illustrated in |
| 1213 |
> |
Fig.~\ref{oopseFig:dynamicPropsMemory}. |
| 1214 |
|
|
| 1215 |
< |
\begin{itemize} |
| 1216 |
< |
|
| 1217 |
< |
\item Restate capabilities |
| 1215 |
> |
\begin{figure} |
| 1216 |
> |
\centering |
| 1217 |
> |
\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 1218 |
> |
\caption[A representation of the block correlations in \texttt{dynamicProps}]{This diagram illustrates the memory management used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.} |
| 1219 |
> |
\label{oopseFig:dynamicPropsMemory} |
| 1220 |
> |
\end{figure} |
| 1221 |
|
|
| 1222 |
< |
\item recap major structure / design choices |
| 1222 |
> |
\subsection{\label{openSource}Open Source and Distribution License} |
| 1223 |
|
|
| 1224 |
< |
\begin{itemize} |
| 1138 |
< |
|
| 1139 |
< |
\item parallel |
| 1140 |
< |
\item symplectic integration |
| 1141 |
< |
\item languages |
| 1224 |
> |
\section{\label{oopseSec:conclusion}Conclusion} |
| 1225 |
|
|
| 1226 |
< |
\end{itemize} |
| 1226 |
> |
We have presented the design and implementation of our open source |
| 1227 |
> |
simulation package {\sc oopse}. The package offers novel |
| 1228 |
> |
capabilities to the field of Molecular Dynamics simulation packages in |
| 1229 |
> |
the form of dipolar force fields, and symplectic integration of rigid |
| 1230 |
> |
body dynamics. It is capable of scaling across multiple processors |
| 1231 |
> |
through the use of MPI. It also implements several integration |
| 1232 |
> |
ensembles allowing the end user control over temperature and |
| 1233 |
> |
pressure. In addition, it is capable of integrating constrained |
| 1234 |
> |
dynamics through both the {\sc rattle} algorithm and the z-constraint |
| 1235 |
> |
method. |
| 1236 |
|
|
| 1237 |
< |
\item How well does it meet the primary goal |
| 1237 |
> |
These features are all brought together in a single open-source |
| 1238 |
> |
development package. This allows researchers to not only benefit from |
| 1239 |
> |
{\sc oopse}, but also contribute to {\sc oopse}'s development as |
| 1240 |
> |
well.Documentation and source code for {\sc oopse} can be downloaded |
| 1241 |
> |
from \texttt{http://www.openscience.org/oopse/}. |
| 1242 |
|
|
| 1147 |
– |
\end{itemize} |
