| 3 |
|
\usepackage{amssymb} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
|
\usepackage{listings} |
| 6 |
< |
\usepackage{palatino} |
| 6 |
> |
\usepackage{berkeley} |
| 7 |
|
\usepackage{graphicx} |
| 8 |
|
\usepackage[ref]{overcite} |
| 9 |
|
\usepackage{setspace} |
| 35 |
|
\maketitle |
| 36 |
|
|
| 37 |
|
\begin{abstract} |
| 38 |
< |
Need an abstract. |
| 38 |
> |
{\sc oopse} is a new molecular dynamics simulation program which is |
| 39 |
> |
capable of efficiently integrating equations of motion for atom types |
| 40 |
> |
with orientational degrees of freedom (e.g. ``sticky'' atoms and point |
| 41 |
> |
dipoles). Transition metals can also be simulated using the embedded |
| 42 |
> |
atom method ({\sc eam}) potential included in the code. Parallel |
| 43 |
> |
simulations are carried out using the force-based decomposition |
| 44 |
> |
method. Simulations are specified using a very simple C-based |
| 45 |
> |
meta-data language. A number of advanced integrators are included, |
| 46 |
> |
and the base integrator for orientational dynamics provides |
| 47 |
> |
substantial improvements over older quaternion-base schemes. All |
| 48 |
> |
source code is available under a very permissive (BSD-style) Open |
| 49 |
> |
Source license. |
| 50 |
|
\end{abstract} |
| 51 |
|
|
| 52 |
|
\section{\label{sec:intro}Introduction} |
| 53 |
|
|
| 54 |
< |
UNDERWAY |
| 54 |
> |
There are a number of excellent molecular dynamics packages available |
| 55 |
> |
to the chemical physics |
| 56 |
> |
community.\cite{Brooks83,MacKerell98,pearlman:1995,Gromacs,Gromacs3,DL_POLY,Tinker,Paradyn} |
| 57 |
> |
All of these packages are stable, polished programs which solve many |
| 58 |
> |
problems of interest. Most are now capable of performing molecular |
| 59 |
> |
dynamics simulations on parallel computers. Some have source code |
| 60 |
> |
which is freely available to the entire scientific community. Few, |
| 61 |
> |
however, are capable of efficiently integrating the equations of |
| 62 |
> |
motion for atom types with orientational degrees of freedom |
| 63 |
> |
(e.g. point dipoles, and ``sticky'' atoms). And only one of the |
| 64 |
> |
programs referenced can handle transition metal force fields like the |
| 65 |
> |
Embedded Atom Method ({\sc eam}). The direction our research program |
| 66 |
> |
has taken us now involves the use of atoms with orientational degrees |
| 67 |
> |
of freedom and transition metals. Since these simulation methods may |
| 68 |
> |
be of some use to other researchers, we have decided to release our |
| 69 |
> |
program to the scientific community with a permissive open source |
| 70 |
> |
license. |
| 71 |
|
|
| 72 |
< |
|
| 73 |
< |
We have structured this paper to first discuss the underlying concepts |
| 74 |
< |
in this simulation package (Sec. \ref{oopseSec:IOfiles}). The |
| 75 |
< |
empirical energy functions implemented are discussed in |
| 76 |
< |
Sec.~\ref{oopseSec:empiricalEnergy}. Sec.~\ref{oopseSec:mechanics} |
| 77 |
< |
describes the various Molecular Dynamics algorithms {\sc oopse} |
| 78 |
< |
implements in the integration of the Newtonian equations of motion. |
| 79 |
< |
Program design considerations for parallel computing are presented in |
| 72 |
> |
This paper communicates the algorithmic details of our program, which |
| 73 |
> |
we have been calling the Open source Object-oriented Parallel |
| 74 |
> |
Simulation Engine (i.e. {\sc oopse}). We have structured this paper |
| 75 |
> |
to first discuss the underlying concepts in this simulation package |
| 76 |
> |
(Sec. \ref{oopseSec:IOfiles}). The empirical energy functions |
| 77 |
> |
implemented are discussed in Sec.~\ref{oopseSec:empiricalEnergy}. |
| 78 |
> |
Sec.~\ref{oopseSec:mechanics} describes the various Molecular Dynamics |
| 79 |
> |
algorithms {\sc oopse} implements in the integration of Hamilton's |
| 80 |
> |
equations of motion. Program design considerations for parallel |
| 81 |
> |
computing are presented in |
| 82 |
|
Sec.~\ref{oopseSec:parallelization}. Concluding remarks are presented |
| 83 |
|
in Sec.~\ref{oopseSec:conclusion}. |
| 84 |
|
|
| 366 |
|
Sec.~\ref{oopseSec:minimizer}. The {\tt forceField} statement is |
| 367 |
|
important for the selection of which forces will be used in the course |
| 368 |
|
of the simulation. {\sc oopse} supports several force fields, as |
| 369 |
< |
outlined in Sec.~\ref{oopseSec:empericalEnergy}. The force fields are |
| 369 |
> |
outlined in Sec.~\ref{oopseSec:empiricalEnergy}. The force fields are |
| 370 |
|
interchangeable between simulations, with the only requirement being |
| 371 |
|
that all atoms needed by the simulation are defined within the |
| 372 |
|
selected force field. |
| 1077 |
|
and the drawbacks and benefits of the different density corrected SSD |
| 1078 |
|
models can be found in reference~\citen{fennell04}. |
| 1079 |
|
|
| 1080 |
+ |
\subsection{\label{oopseSec:WATER}The {\sc water} Force Field} |
| 1081 |
+ |
|
| 1082 |
+ |
In addition to the {\sc duff} force field's solvent description, a |
| 1083 |
+ |
separate {\sc water} force field has been included for simulating many |
| 1084 |
+ |
of the common rigid-body water models. In addition to the simple or |
| 1085 |
+ |
dipolar models (SSD, SSD1, SSD/E, SSD/RF, and DPD water), the common |
| 1086 |
+ |
charge-based models were included (SPC, SPC/E, TIP3P, TIP4P, and |
| 1087 |
+ |
TIP5P).\cite{liu96:new_model,Ichiye03,fennell04,Marrink01,Berendsen81,Berendsen87,Jorgensen83,Mahoney00} |
| 1088 |
+ |
In order to handle these models, charge-charge interactions were |
| 1089 |
+ |
included in the force-loop: |
| 1090 |
+ |
\begin{equation} |
| 1091 |
+ |
V_{\text{charge}}(r_{ij}) = \sum_{ij}\frac{q_iq_je^2}{r_{ij}}, |
| 1092 |
+ |
\end{equation} |
| 1093 |
+ |
where $q$ represents the charge on particle $i$ or $j$, and $e$ is the |
| 1094 |
+ |
charge of an electron in Coulombs. The charge-charge interaction |
| 1095 |
+ |
support is rudimentary in the current version of {\sc oopse}. As with |
| 1096 |
+ |
the other pair interactions, charges can be simulated with a pure |
| 1097 |
+ |
cutoff or a reaction field. The various methods for performing the |
| 1098 |
+ |
Ewald summation have not yet been included. Also, the charge-dipole |
| 1099 |
+ |
and charge-quadrupole (for interactions between SSD type water and |
| 1100 |
+ |
charges) are not yet available, so it is currently inadvisable to mix |
| 1101 |
+ |
dipolar and charge based molecules in the same system. |
| 1102 |
+ |
|
| 1103 |
+ |
The {\sc water} force field can be easily expanded through |
| 1104 |
+ |
modification of the {\sc water} force field file ({\tt WATER.frc}). By |
| 1105 |
+ |
adding atom types and inserting the appropriate parameters, it is |
| 1106 |
+ |
possible to extend the force field to handle rigid molecules other |
| 1107 |
+ |
than water. |
| 1108 |
+ |
|
| 1109 |
|
\subsection{\label{oopseSec:eam}Embedded Atom Method} |
| 1110 |
|
|
| 1111 |
|
{\sc oopse} implements a potential that describes bonding in |
| 1467 |
|
|
| 1468 |
|
The matrix rotations used in the DLM method end up being more costly |
| 1469 |
|
computationally than the simpler arithmetic quaternion |
| 1470 |
< |
propagation. With the same time step, a 1000-molecule water simulation |
| 1471 |
< |
shows an average 7\% increase in computation time using the DLM method |
| 1472 |
< |
in place of quaternions. This cost is more than justified when |
| 1473 |
< |
comparing the energy conservation of the two methods as illustrated in |
| 1474 |
< |
Fig.~\ref{timestep}. |
| 1470 |
> |
propagation. With the same time step, a 1024-molecule water simulation |
| 1471 |
> |
shows an 12\% increase in computation time (averaged over several |
| 1472 |
> |
different time steps) using the DLM method in place of |
| 1473 |
> |
quaternions. This cost is more than justified when comparing the |
| 1474 |
> |
energy conservation of the two methods. Figure ~\ref{quatdlm} provides |
| 1475 |
> |
a comparative analysis of the {\sc dlm} method versus the simple |
| 1476 |
> |
quaternion method that was originally implemented. |
| 1477 |
|
|
| 1478 |
|
\begin{figure} |
| 1479 |
|
\centering |
| 1480 |
< |
\includegraphics[width=\linewidth]{timeStep.eps} |
| 1481 |
< |
\caption[Energy conservation for quaternion versus DLM dynamics]{Energy conservation using quaternion based integration versus |
| 1482 |
< |
the method proposed by Dullweber \emph{et al.} with increasing time |
| 1483 |
< |
step. For each time step, the dotted line is total energy using the |
| 1484 |
< |
DLM integrator, and the solid line comes from the quaternion |
| 1485 |
< |
integrator. The larger time step plots are shifted up from the true |
| 1486 |
< |
energy baseline for clarity.} |
| 1487 |
< |
\label{timestep} |
| 1480 |
> |
\includegraphics[width=\linewidth]{quatvsdlm.eps} |
| 1481 |
> |
\caption[Energy conservation analysis of the {\sc dlm} and quaternion |
| 1482 |
> |
integration methods]{The logarithm of absolute value of the slope of |
| 1483 |
> |
the energy drift (\delta E$_1$) and the standard deviation of the |
| 1484 |
> |
energy fluctuations (\delta E$_0$) as a function of chosen time |
| 1485 |
> |
step. All simulations were of a 1024-molecule simulation of SSD water |
| 1486 |
> |
at 298 K starting from the same initial configuration. Note that the |
| 1487 |
> |
{\sc dlm} method provides a greater-than order-of-magnitude |
| 1488 |
> |
improvement in energy conservation and relative energy fluctuations |
| 1489 |
> |
over the quaternion method at all the tested time steps. The energy |
| 1490 |
> |
drift is quite steep for the larger time steps in both methods, and |
| 1491 |
> |
results in discontinuous behavior as the systems compound their |
| 1492 |
> |
anomolous energy accumulation.} |
| 1493 |
> |
\label{quatdlm} |
| 1494 |
|
\end{figure} |
| 1495 |
|
|
| 1496 |
< |
In Fig.~\ref{timestep}, the resulting energy drift at various time |
| 1497 |
< |
steps for both the DLM and quaternion integration schemes is |
| 1498 |
< |
compared. All of the 1000 molecule water simulations started with the |
| 1499 |
< |
same configuration, and the only difference was the method for |
| 1500 |
< |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 1501 |
< |
methods for propagating molecule rotation conserve energy fairly well, |
| 1502 |
< |
with the quaternion method showing a slight energy drift over time in |
| 1503 |
< |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 1504 |
< |
energy conservation benefits of the DLM method are clearly |
| 1505 |
< |
demonstrated. Thus, while maintaining the same degree of energy |
| 1506 |
< |
conservation, one can take considerably longer time steps, leading to |
| 1507 |
< |
an overall reduction in computation time. |
| 1496 |
> |
In Fig.~\ref{quatdlm}, \delta E$_1$ is a measure of the linear energy |
| 1497 |
> |
drift in units of kcal/mol per particle over a nanosecond of |
| 1498 |
> |
simulation time, and \delta E$_0$ is the standard deviation of the |
| 1499 |
> |
energy fluctuations in units of kcal/mol per particle. In the top |
| 1500 |
> |
plot, it is apparent that the energy drift is reduced by a significant |
| 1501 |
> |
amount (2 to 3 orders-of-magnitude improvement at every tested time |
| 1502 |
> |
step) by chosing the {\sc dlm} method over the simple non-symplectic |
| 1503 |
> |
quaternion integration method. When the energy drift becomes very |
| 1504 |
> |
small ($log_{10}[|\delta\text{E}_1|] < -3$), it is more difficult to |
| 1505 |
> |
calculate a slope, resulting in the larger displayed error bars. In |
| 1506 |
> |
addition to this improvement in energy drift, the fluctuation is the |
| 1507 |
> |
total energy are also dampened out by 1 to 2 orders-of-magnitude by |
| 1508 |
> |
utilizing the {\sc dlm} integration method. |
| 1509 |
|
|
| 1510 |
+ |
It was stated previously that the {\sc dlm} method was the more |
| 1511 |
+ |
computationally expensive of the two implimented integration |
| 1512 |
+ |
methodologies. In order to incorporate this information into the |
| 1513 |
+ |
energy analysis a plot of energy drift versus computational cost was |
| 1514 |
+ |
generated (Fig.~\ref{cpuCost}). This figure provides an estimate of |
| 1515 |
+ |
the CPU time required under the two integration schemes for 1 |
| 1516 |
+ |
nanosecond of simulation time for the model 1024-molecule system. The |
| 1517 |
+ |
plot is read by chosing a desired energy drift value and determining |
| 1518 |
+ |
where both the curves cross. If a \delta E$_1$ of 1E-3 kcal/mol per |
| 1519 |
+ |
particle is desired, a nanosecond of simulation time will require ~19 |
| 1520 |
+ |
hours of CPU time with the {\sc dlm} integrator, while the same small |
| 1521 |
+ |
drift value will require ~154 hours of CPU time. This demonstrates the |
| 1522 |
+ |
computational advantage of the integration scheme utilized in {\sc |
| 1523 |
+ |
oopse}. |
| 1524 |
+ |
|
| 1525 |
+ |
\begin{figure} |
| 1526 |
+ |
\centering |
| 1527 |
+ |
\includegraphics[width=\linewidth]{compCost.eps} |
| 1528 |
+ |
\caption[Energy drift as a function of required simulation run |
| 1529 |
+ |
time]{The logarithm of absolute value of the slope of the energy drift |
| 1530 |
+ |
(\delta E$_1$) as a function of simulation run time. Simulations were |
| 1531 |
+ |
performed on a single 2.5 GHz Pentium IV processor. Simulation time |
| 1532 |
+ |
comparisons can be made by tracing horizontally from one curve to the |
| 1533 |
+ |
other. For example, a simulation that takes ~24 hours using the {\sc |
| 1534 |
+ |
dlm} method will take roughly 210 hours using the simple quaternion |
| 1535 |
+ |
method if the same degree of energy conservation is desired.} |
| 1536 |
+ |
\label{cpuCost} |
| 1537 |
+ |
\end{figure} |
| 1538 |
+ |
|
| 1539 |
|
There is only one specific keyword relevant to the default integrator, |
| 1540 |
|
and that is the time step for integrating the equations of motion. |
| 1541 |
|
|
| 2177 |
|
\end{table} |
| 2178 |
|
|
| 2179 |
|
|
| 2180 |
< |
\section{\label{sec:minimize}Energy Minimization} |
| 2180 |
> |
\section{\label{oopseSec:minimizer}Energy Minimization} |
| 2181 |
|
|
| 2182 |
|
As one of the basic procedures of molecular modeling, energy |
| 2183 |
|
minimization is used to identify local configurations that are stable |
