| 2 |
|
\usepackage{amsmath} |
| 3 |
|
\usepackage{amssymb} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
– |
%\usepackage{berkeley} |
| 5 |
|
\usepackage{listings} |
| 6 |
+ |
\usepackage{palatino} |
| 7 |
|
\usepackage{graphicx} |
| 8 |
|
\usepackage[ref]{overcite} |
| 9 |
|
\usepackage{setspace} |
| 77 |
|
those of the community of scientists who contribute back to the |
| 78 |
|
project. |
| 79 |
|
|
| 80 |
< |
We have structured this chapter to first discuss the empirical energy |
| 80 |
> |
We have structured this paper to first discuss the empirical energy |
| 81 |
|
functions that {\sc oopse } implements in |
| 82 |
|
Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of |
| 83 |
|
the various input and output files associated with the package |
| 84 |
|
(Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics} |
| 85 |
|
elucidates the various Molecular Dynamics algorithms {\sc oopse} |
| 86 |
|
implements in the integration of the Newtonian equations of |
| 87 |
< |
motion. Basic analysis of the trajectories obtained from the |
| 88 |
< |
simulation is discussed in Sec.~\ref{oopseSec:props}. Program design |
| 87 |
> |
motion. Program design |
| 88 |
|
considerations are presented in Sec.~\ref{oopseSec:design}. And |
| 89 |
|
lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter. |
| 90 |
|
|
| 795 |
|
entities are written out using quanternions, to save space in the |
| 796 |
|
output files. |
| 797 |
|
|
| 798 |
< |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular velocities.},label=sch:dumpFormat] |
| 798 |
> |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly, body fixed angular momentum.},label=sch:dumpFormat] |
| 799 |
|
|
| 800 |
|
nAtoms |
| 801 |
|
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 809 |
|
There are three major files used by {\sc oopse} written in the |
| 810 |
|
coordinate format, they are as follows: the initialization file |
| 811 |
|
(\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and |
| 812 |
< |
the final coordinates of the simulation. The initialization file is |
| 812 |
> |
the final coordinates of the simulation (\texttt{.eor}). The initialization file is |
| 813 |
|
necessary for {\sc oopse} to start the simulation with the proper |
| 814 |
|
coordinates, and is generated before the simulation run. The |
| 815 |
|
trajectory file is created at the beginning of the simulation, and is |
| 830 |
|
|
| 831 |
|
As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization |
| 832 |
|
file is needed to provide the starting coordinates for a |
| 833 |
< |
simulation. The {\sc oopse} package provides several system building |
| 834 |
< |
programs to aid in the creation of the \texttt{.init} |
| 835 |
< |
file. The programs use {\sc bass}, and will recognize |
| 836 |
< |
arguments and parameters in the \texttt{.bass} file that would |
| 837 |
< |
otherwise be ignored by the simulation. |
| 833 |
> |
simulation. Several helper programs are provided with {\sc oopse} to illustrate possible build routes. However, as each simulation is different, system creation is left to the end user. The {\tt .init} file must list the atoms in the correct order or {\sc oopse} will give an atom mismatch error. |
| 834 |
> |
|
| 835 |
> |
The correct ordering of the atoms relies on the ordering of atoms and molecules within the model and {\sc bass} scripts. {\sc oopse} expects the order to comply with the following guidelines: |
| 836 |
> |
\begin{enumerate} |
| 837 |
> |
\item All of the molecules of the first declared component are given before proceeding to the molecules of the second component, and so on for all declared components. |
| 838 |
> |
\item The ordering of the atoms for each molecule follows the order declared in the molecule's declaration within the model file. |
| 839 |
> |
\end{enumerate} |
| 840 |
> |
An example is given in Scheme~\ref{sch:initEx1} resulting in the {\tt .init} file shown in Scheme~\ref{sch:initEx2}. |
| 841 |
> |
|
| 842 |
> |
\begin{lstlisting}[float,caption={This scheme illustrates the declaration of the $\text{I}_2$ molecule and the HCl molecule. The two molecules are then included into a simulation.}, label=sch:initEx1] |
| 843 |
> |
|
| 844 |
> |
molecule{ |
| 845 |
> |
name = "I2"; |
| 846 |
> |
nAtoms = 2; |
| 847 |
> |
atom[0]{ |
| 848 |
> |
type = "I"; |
| 849 |
> |
} |
| 850 |
> |
atom[1]{ |
| 851 |
> |
type = "I"; |
| 852 |
> |
} |
| 853 |
> |
nBonds = 1; |
| 854 |
> |
bond[0]{ |
| 855 |
> |
members( 0, 1); |
| 856 |
> |
} |
| 857 |
> |
} |
| 858 |
|
|
| 859 |
+ |
molecule{ |
| 860 |
+ |
name = "HCl" |
| 861 |
+ |
nAtoms = 2; |
| 862 |
+ |
atom[0]{ |
| 863 |
+ |
type = "H"; |
| 864 |
+ |
} |
| 865 |
+ |
atom[1]{ |
| 866 |
+ |
type = "Cl"; |
| 867 |
+ |
} |
| 868 |
+ |
nBonds = 1; |
| 869 |
+ |
bond[0]{ |
| 870 |
+ |
members( 0, 1); |
| 871 |
+ |
} |
| 872 |
+ |
} |
| 873 |
+ |
|
| 874 |
+ |
nComponents = 2; |
| 875 |
+ |
component{ |
| 876 |
+ |
type = "HCl"; |
| 877 |
+ |
nMol = 4; |
| 878 |
+ |
} |
| 879 |
+ |
component{ |
| 880 |
+ |
type = "I2"; |
| 881 |
+ |
nMol = 1; |
| 882 |
+ |
} |
| 883 |
+ |
|
| 884 |
+ |
initialConfig = "mixture.init"; |
| 885 |
+ |
|
| 886 |
+ |
\end{lstlisting} |
| 887 |
+ |
|
| 888 |
+ |
\begin{lstlisting}[float,caption={This is the contents of the {\tt mixture.init} file matching the declarations in Scheme~\ref{sch:initEx1}. Note that even though $\text{I}_2$ is declared before HCl, the {\tt .init} file follows the order in which the components were included.},label=sch:initEx2] |
| 889 |
+ |
|
| 890 |
+ |
10 |
| 891 |
+ |
0.0; 10.0 0.0 0.0; 0.0 10.0 0.0; 0.0 0.0 10.0; |
| 892 |
+ |
H ... |
| 893 |
+ |
Cl ... |
| 894 |
+ |
H ... |
| 895 |
+ |
Cl ... |
| 896 |
+ |
H ... |
| 897 |
+ |
Cl ... |
| 898 |
+ |
H ... |
| 899 |
+ |
Cl ... |
| 900 |
+ |
I ... |
| 901 |
+ |
I ... |
| 902 |
+ |
|
| 903 |
+ |
\end{lstlisting} |
| 904 |
+ |
|
| 905 |
+ |
|
| 906 |
|
\subsection{The Statistics File} |
| 907 |
|
|
| 908 |
|
The last output file generated by {\sc oopse} is the statistics |
| 1703 |
|
orientational anisotropy in the system (i.e. in lipid bilayer |
| 1704 |
|
simulations). |
| 1705 |
|
|
| 1706 |
< |
\subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 1706 |
> |
\subsection{\label{sec:constraints}Constraint Methods} |
| 1707 |
> |
|
| 1708 |
> |
\subsubsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond |
| 1709 |
|
Constraints} |
| 1710 |
|
|
| 1711 |
|
In order to satisfy the constraints of fixed bond lengths within {\sc |
| 1712 |
|
oopse}, we have implemented the {\sc rattle} algorithm of |
| 1713 |
|
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
| 1714 |
|
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 1715 |
< |
solving the Lagrange multipliers of constraint. The system of Lagrange |
| 1648 |
< |
multipliers allows one to reformulate the equations of motion with |
| 1649 |
< |
explicit constraint forces.\cite{fowles99:lagrange} |
| 1650 |
< |
|
| 1651 |
< |
Consider a system described by coordinates $q_1$ and $q_2$ subject to an |
| 1652 |
< |
equation of constraint: |
| 1653 |
< |
\begin{equation} |
| 1654 |
< |
\sigma(q_1, q_2,t) = 0 |
| 1655 |
< |
\label{oopseEq:lm1} |
| 1656 |
< |
\end{equation} |
| 1657 |
< |
The Lagrange formulation of the equations of motion can be written: |
| 1658 |
< |
\begin{equation} |
| 1659 |
< |
\delta\int_{t_1}^{t_2}L\, dt = |
| 1660 |
< |
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
| 1661 |
< |
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
| 1662 |
< |
\biggr ) \biggr] \delta q_i \, dt = 0. |
| 1663 |
< |
\label{oopseEq:lm2} |
| 1664 |
< |
\end{equation} |
| 1665 |
< |
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
| 1666 |
< |
are linked by $\sigma$. However, $\sigma$ is fixed at any given |
| 1667 |
< |
instant of time, giving: |
| 1668 |
< |
\begin{align} |
| 1669 |
< |
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
| 1670 |
< |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 ,\\ |
| 1671 |
< |
% |
| 1672 |
< |
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
| 1673 |
< |
- \frac{\partial\sigma}{\partial q_2} \delta q_2, \\ |
| 1674 |
< |
% |
| 1675 |
< |
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1676 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1. |
| 1677 |
< |
\end{align} |
| 1678 |
< |
Substituted back into Eq.~\ref{oopseEq:lm2}, |
| 1679 |
< |
\begin{equation} |
| 1680 |
< |
\int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1} |
| 1681 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1682 |
< |
\biggr) |
| 1683 |
< |
- \biggl( \frac{\partial L}{\partial q_1} |
| 1684 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1685 |
< |
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1686 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0. |
| 1687 |
< |
\label{oopseEq:lm3} |
| 1688 |
< |
\end{equation} |
| 1689 |
< |
Leading to, |
| 1690 |
< |
\begin{equation} |
| 1691 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1692 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1693 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
| 1694 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_2} |
| 1695 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
| 1696 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_2}}. |
| 1697 |
< |
\label{oopseEq:lm4} |
| 1698 |
< |
\end{equation} |
| 1699 |
< |
This relation can only be statisfied, if both are equal to a single |
| 1700 |
< |
function $-\lambda(t)$, |
| 1701 |
< |
\begin{align} |
| 1702 |
< |
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1703 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1704 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t), \\ |
| 1705 |
< |
% |
| 1706 |
< |
\frac{\partial L}{\partial q_1} |
| 1707 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
| 1708 |
< |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} ,\\ |
| 1709 |
< |
% |
| 1710 |
< |
\frac{\partial L}{\partial q_1} |
| 1711 |
< |
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1712 |
< |
+ \mathcal{G}_i &= 0, |
| 1713 |
< |
\end{align} |
| 1714 |
< |
where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1715 |
< |
\begin{equation} |
| 1716 |
< |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}. |
| 1717 |
< |
\label{oopseEq:lm5} |
| 1718 |
< |
\end{equation} |
| 1715 |
> |
solving the Lagrange multipliers of constraint. |
| 1716 |
|
|
| 1717 |
< |
In a simulation, this would involve the solution of a set of $(m + n)$ |
| 1721 |
< |
number of equations. Where $m$ is the number of constraints, and $n$ |
| 1722 |
< |
is the number of constrained coordinates. In practice, this is not |
| 1723 |
< |
done, as the matrix inversion necessary to solve the system of |
| 1724 |
< |
equations would be very time consuming to solve. Additionally, the |
| 1725 |
< |
numerical error in the solution of the set of $\lambda$'s would be |
| 1726 |
< |
compounded by the error inherent in propagating by the Velocity Verlet |
| 1727 |
< |
algorithm ($\Delta t^4$). The Verlet propagation error is negligible |
| 1728 |
< |
in an unconstrained system, as one is interested in the statistics of |
| 1729 |
< |
the run, and not that the run be numerically exact to the ``true'' |
| 1730 |
< |
integration. This relates back to the ergodic hypothesis that a time |
| 1731 |
< |
integral of a valid trajectory will still give the correct ensemble |
| 1732 |
< |
average. However, in the case of constraints, if the equations of |
| 1733 |
< |
motion leave the ``true'' trajectory, they are departing from the |
| 1734 |
< |
constrained surface. The method that is used, is to iteratively solve |
| 1735 |
< |
for $\lambda(t)$ at each time step. |
| 1736 |
< |
|
| 1737 |
< |
In {\sc rattle} the equations of motion are modified subject to the |
| 1738 |
< |
following two constraints: |
| 1739 |
< |
\begin{align} |
| 1740 |
< |
\sigma_{ij}[\mathbf{r}(t)] \equiv |
| 1741 |
< |
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
| 1742 |
< |
\label{oopseEq:c1}, \\ |
| 1743 |
< |
% |
| 1744 |
< |
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
| 1745 |
< |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 .\label{oopseEq:c2} |
| 1746 |
< |
\end{align} |
| 1747 |
< |
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1748 |
< |
the constrained distance between atom $i$ and |
| 1749 |
< |
$j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to |
| 1750 |
< |
be perpendicular to the bond vector, so that the bond can neither grow |
| 1751 |
< |
nor shrink. The constrained dynamics equations become: |
| 1752 |
< |
\begin{equation} |
| 1753 |
< |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i, |
| 1754 |
< |
\label{oopseEq:r1} |
| 1755 |
< |
\end{equation} |
| 1756 |
< |
where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$, |
| 1757 |
< |
and are defined: |
| 1758 |
< |
\begin{equation} |
| 1759 |
< |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}. |
| 1760 |
< |
\label{oopseEq:r2} |
| 1761 |
< |
\end{equation} |
| 1762 |
< |
|
| 1763 |
< |
In Velocity Verlet, if $\Delta t = h$, the propagation can be written: |
| 1764 |
< |
\begin{align} |
| 1765 |
< |
\mathbf{r}_i(t+h) &= |
| 1766 |
< |
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
| 1767 |
< |
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
| 1768 |
< |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1}, \\ |
| 1769 |
< |
% |
| 1770 |
< |
\mathbf{\dot{r}}_i(t+h) &= |
| 1771 |
< |
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
| 1772 |
< |
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
| 1773 |
< |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] ,% |
| 1774 |
< |
\label{oopseEq:vv2} |
| 1775 |
< |
\end{align} |
| 1776 |
< |
where: |
| 1777 |
< |
\begin{align} |
| 1778 |
< |
\mathbf{\mathcal{G}}_{Ri}(t) &= |
| 1779 |
< |
-2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) ,\\ |
| 1780 |
< |
% |
| 1781 |
< |
\mathbf{\mathcal{G}}_{Vi}(t+h) &= |
| 1782 |
< |
-2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h). |
| 1783 |
< |
\end{align} |
| 1784 |
< |
Next, define: |
| 1785 |
< |
\begin{align} |
| 1786 |
< |
g_{ij} &= h \lambda_{Rij}(t) ,\\ |
| 1787 |
< |
k_{ij} &= h \lambda_{Vij}(t+h), \\ |
| 1788 |
< |
\mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t) |
| 1789 |
< |
- \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t). |
| 1790 |
< |
\end{align} |
| 1791 |
< |
Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2} |
| 1792 |
< |
can be rewritten as, |
| 1793 |
< |
\begin{align} |
| 1794 |
< |
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i ,\\ |
| 1795 |
< |
% |
| 1796 |
< |
\mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1797 |
< |
-\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h). |
| 1798 |
< |
\end{align} |
| 1799 |
< |
|
| 1800 |
< |
To integrate the equations of motion, the {\sc rattle} algorithm first |
| 1801 |
< |
solves for $\mathbf{r}(t+h)$. Let, |
| 1802 |
< |
\begin{equation} |
| 1803 |
< |
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t). |
| 1804 |
< |
\end{equation} |
| 1805 |
< |
Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next |
| 1806 |
< |
pick a constraint $j$, and let, |
| 1807 |
< |
\begin{equation} |
| 1808 |
< |
\mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t) |
| 1809 |
< |
- \mathbf{r}_j(t) + h\mathbf{q}_j(t). |
| 1810 |
< |
\label{oopseEq:ra1} |
| 1811 |
< |
\end{equation} |
| 1812 |
< |
If |
| 1813 |
< |
\begin{equation} |
| 1814 |
< |
\Big| |\mathbf{s}|^2 - d_{ij}^2 \Big| > \text{tolerance}, |
| 1815 |
< |
\end{equation} |
| 1816 |
< |
then the constraint is unsatisfied, and corrections are made to the |
| 1817 |
< |
positions. First we define a test corrected configuration as, |
| 1818 |
< |
\begin{align} |
| 1819 |
< |
\mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i - |
| 1820 |
< |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] ,\\ |
| 1821 |
< |
% |
| 1822 |
< |
\mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j + |
| 1823 |
< |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr]. |
| 1824 |
< |
\end{align} |
| 1825 |
< |
And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2 |
| 1826 |
< |
= d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the |
| 1827 |
< |
approximation, |
| 1828 |
< |
\begin{equation} |
| 1829 |
< |
g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)] |
| 1830 |
< |
(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1831 |
< |
\end{equation} |
| 1832 |
< |
Although not an exact solution for $g_{ij}$, as this is an iterative |
| 1833 |
< |
scheme overall, the eventual solution will converge. With a trial |
| 1834 |
< |
$g_{ij}$, the new $\mathbf{q}$'s become, |
| 1835 |
< |
\begin{align} |
| 1836 |
< |
\mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\, |
| 1837 |
< |
\frac{\mathbf{r}_{ij}(t)}{m_i} ,\\ |
| 1838 |
< |
% |
| 1839 |
< |
\mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\, |
| 1840 |
< |
\frac{\mathbf{r}_{ij}(t)}{m_j} . |
| 1841 |
< |
\end{align} |
| 1842 |
< |
The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until |
| 1843 |
< |
all constraints are satisfied. |
| 1717 |
> |
\subsubsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1718 |
|
|
| 1845 |
– |
The second step of {\sc rattle}, is to then update the velocities. The |
| 1846 |
– |
step starts with, |
| 1847 |
– |
\begin{equation} |
| 1848 |
– |
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h). |
| 1849 |
– |
\end{equation} |
| 1850 |
– |
Next we pick a constraint $j$, and calculate the dot product $\ell$. |
| 1851 |
– |
\begin{equation} |
| 1852 |
– |
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h). |
| 1853 |
– |
\label{oopseEq:rv1} |
| 1854 |
– |
\end{equation} |
| 1855 |
– |
Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be |
| 1856 |
– |
zero. Therefore if $\ell$ is greater than some tolerance, then |
| 1857 |
– |
corrections are made to the $i$ and $j$ velocities. |
| 1858 |
– |
\begin{align} |
| 1859 |
– |
\mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij} |
| 1860 |
– |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i}, \\ |
| 1861 |
– |
% |
| 1862 |
– |
\mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij} |
| 1863 |
– |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j}. |
| 1864 |
– |
\end{align} |
| 1865 |
– |
Like in the previous step, we select a value for $k_{ij}$ such that |
| 1866 |
– |
$\ell$ is zero. |
| 1867 |
– |
\begin{equation} |
| 1868 |
– |
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1869 |
– |
\end{equation} |
| 1870 |
– |
The test velocities, $\mathbf{\dot{r}}^T_i$ and |
| 1871 |
– |
$\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and |
| 1872 |
– |
the algorithm is iterated from Eq.~\ref{oopseEq:rv1} until all |
| 1873 |
– |
constraints are satisfied. |
| 1874 |
– |
|
| 1875 |
– |
|
| 1876 |
– |
\subsection{\label{oopseSec:zcons}Z-Constraint Method} |
| 1877 |
– |
|
| 1719 |
|
Based on the fluctuation-dissipation theorem, a force auto-correlation |
| 1720 |
|
method was developed by Roux and Karplus to investigate the dynamics |
| 1721 |
|
of ions inside ion channels.\cite{Roux91} The time-dependent friction |
| 1802 |
|
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 1803 |
|
\end{equation} |
| 1804 |
|
|
| 1964 |
– |
\section{\label{oopseSec:props}Trajectory Analysis} |
| 1965 |
– |
|
| 1966 |
– |
\subsection{\label{oopseSec:staticProps}Static Property Analysis} |
| 1967 |
– |
|
| 1968 |
– |
The static properties of the trajectories are analyzed with the |
| 1969 |
– |
program \texttt{staticProps}. The code is capable of calculating a |
| 1970 |
– |
number of pair correlations between species A and B. Some of which |
| 1971 |
– |
only apply to directional entities. The summary of pair correlations |
| 1972 |
– |
can be found in Table~\ref{oopseTb:gofrs} |
| 1973 |
– |
|
| 1974 |
– |
\begin{table} |
| 1975 |
– |
\caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}} |
| 1976 |
– |
\label{oopseTb:gofrs} |
| 1977 |
– |
\begin{center} |
| 1978 |
– |
\begin{tabular}{|l|c|c|} |
| 1979 |
– |
\hline |
| 1980 |
– |
Name & Equation & Directional Atom \\ \hline |
| 1981 |
– |
$g_{\text{AB}}(r)$ & Eq.~\ref{eq:gofr} & neither \\ \hline |
| 1982 |
– |
$g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline |
| 1983 |
– |
$g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline |
| 1984 |
– |
$g_{\text{AB}}(x, y, z)$ & Eq.~\ref{eq:gofrXYZ} & neither \\ \hline |
| 1985 |
– |
$\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &% |
| 1986 |
– |
both \\ \hline |
| 1987 |
– |
\end{tabular} |
| 1988 |
– |
\begin{minipage}{\linewidth} |
| 1989 |
– |
\centering |
| 1990 |
– |
\vspace{2mm} |
| 1991 |
– |
The third column specifies which atom, if any, need be a directional entity. |
| 1992 |
– |
\end{minipage} |
| 1993 |
– |
\end{center} |
| 1994 |
– |
\end{table} |
| 1995 |
– |
|
| 1996 |
– |
The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows: |
| 1997 |
– |
\begin{equation} |
| 1998 |
– |
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %% |
| 1999 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} %% |
| 2000 |
– |
\delta( r - |\mathbf{r}_{ij}|) \rangle, \label{eq:gofr} |
| 2001 |
– |
\end{equation} |
| 2002 |
– |
where $\mathbf{r}_{ij}$ is the vector |
| 2003 |
– |
\begin{equation*} |
| 2004 |
– |
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i, \notag |
| 2005 |
– |
\end{equation*} |
| 2006 |
– |
and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over |
| 2007 |
– |
the expected pair density at a given $r$. |
| 2008 |
– |
|
| 2009 |
– |
The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and |
| 2010 |
– |
$g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two |
| 2011 |
– |
dimensional histograms. Both use $r$ for the primary axis then a |
| 2012 |
– |
$\cos$ for the secondary axis ($\cos \theta$ for |
| 2013 |
– |
Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for |
| 2014 |
– |
Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to |
| 2015 |
– |
correlate alignment on directional entities. $g_{\text{AB}}(r, \cos |
| 2016 |
– |
\theta)$ is defined as follows: |
| 2017 |
– |
\begin{equation} |
| 2018 |
– |
g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2019 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2020 |
– |
\delta( \cos \theta - \cos \theta_{ij}) |
| 2021 |
– |
\delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2022 |
– |
\label{eq:gofrCosTheta} |
| 2023 |
– |
\end{equation} |
| 2024 |
– |
Here |
| 2025 |
– |
\begin{equation*} |
| 2026 |
– |
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij}, |
| 2027 |
– |
\end{equation*} |
| 2028 |
– |
where $\mathbf{\hat{i}}$ is the unit directional vector of species $i$ |
| 2029 |
– |
and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector |
| 2030 |
– |
$\mathbf{r}_{ij}$. |
| 2031 |
– |
|
| 2032 |
– |
The second two dimensional histogram is of the form: |
| 2033 |
– |
\begin{equation} |
| 2034 |
– |
g_{\text{AB}}(r, \cos \omega) = |
| 2035 |
– |
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2036 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2037 |
– |
\delta( \cos \omega - \cos \omega_{ij}) |
| 2038 |
– |
\delta( r - |\mathbf{r}_{ij}|) \rangle. \label{eq:gofrCosOmega} |
| 2039 |
– |
\end{equation} |
| 2040 |
– |
Here |
| 2041 |
– |
\begin{equation*} |
| 2042 |
– |
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}}. |
| 2043 |
– |
\end{equation*} |
| 2044 |
– |
Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit |
| 2045 |
– |
directional vectors of species $i$ and $j$. |
| 2046 |
– |
|
| 2047 |
– |
The static analysis code is also cable of calculating a three |
| 2048 |
– |
dimensional pair correlation of the form: |
| 2049 |
– |
\begin{equation}\label{eq:gofrXYZ} |
| 2050 |
– |
g_{\text{AB}}(x, y, z) = |
| 2051 |
– |
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2052 |
– |
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2053 |
– |
\delta( x - x_{ij}) |
| 2054 |
– |
\delta( y - y_{ij}) |
| 2055 |
– |
\delta( z - z_{ij}) \rangle, |
| 2056 |
– |
\end{equation} |
| 2057 |
– |
where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$ |
| 2058 |
– |
components respectively of vector $\mathbf{r}_{ij}$. |
| 2059 |
– |
|
| 2060 |
– |
The final pair correlation is similar to |
| 2061 |
– |
Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega |
| 2062 |
– |
\rangle_{\text{AB}}(r)$ is calculated in the following way: |
| 2063 |
– |
\begin{equation}\label{eq:cosOmegaOfR} |
| 2064 |
– |
\langle \cos \omega \rangle_{\text{AB}}(r) = |
| 2065 |
– |
\langle \sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2066 |
– |
(\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2067 |
– |
\end{equation} |
| 2068 |
– |
Here $\cos \omega_{ij}$ is defined in the same way as in |
| 2069 |
– |
Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair |
| 2070 |
– |
correlation that gives the average correlation of two directional |
| 2071 |
– |
entities as a function of their distance from each other. |
| 2072 |
– |
|
| 2073 |
– |
\subsection{\label{dynamicProps}Dynamic Property Analysis} |
| 2074 |
– |
|
| 2075 |
– |
The dynamic properties of a trajectory are calculated with the program |
| 2076 |
– |
\texttt{dynamicProps}. The program calculates the following properties: |
| 2077 |
– |
\begin{gather} |
| 2078 |
– |
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle, \label{eq:rms}\\ |
| 2079 |
– |
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle, \label{eq:velCorr} \\ |
| 2080 |
– |
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle. \label{eq:angularVelCorr} |
| 2081 |
– |
\end{gather} |
| 2082 |
– |
|
| 2083 |
– |
Eq.~\ref{eq:rms} is the root mean square displacement function. Which |
| 2084 |
– |
allows one to observe the average displacement of an atom as a |
| 2085 |
– |
function of time. The quantity is useful when calculating diffusion |
| 2086 |
– |
coefficients because of the Einstein Relation, which is valid at long |
| 2087 |
– |
times.\cite{allen87:csl} |
| 2088 |
– |
\begin{equation} |
| 2089 |
– |
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle. |
| 2090 |
– |
\label{oopseEq:einstein} |
| 2091 |
– |
\end{equation} |
| 2092 |
– |
|
| 2093 |
– |
Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational |
| 2094 |
– |
velocity and angular velocity correlation functions respectively. The |
| 2095 |
– |
latter is only applicable to directional species in the |
| 2096 |
– |
simulation. The velocity autocorrelation functions are useful when |
| 2097 |
– |
determining vibrational information about the system of interest. |
| 2098 |
– |
|
| 1805 |
|
\section{\label{oopseSec:design}Program Design} |
| 1806 |
|
|
| 1807 |
|
\subsection{\label{sec:architecture} {\sc oopse} Architecture} |
| 1899 |
|
and favorably compete with spatial methods up to 100,000 |
| 1900 |
|
atoms.\cite{plimpton95} |
| 1901 |
|
|
| 2196 |
– |
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 2197 |
– |
|
| 2198 |
– |
For large simulations, the trajectory files can sometimes reach sizes |
| 2199 |
– |
in excess of several gigabytes. In order to effectively analyze that |
| 2200 |
– |
amount of data, two memory management schemes have been devised for |
| 2201 |
– |
\texttt{staticProps} and for \texttt{dynamicProps}. The first scheme, |
| 2202 |
– |
developed for \texttt{staticProps}, is the simplest. As each frame's |
| 2203 |
– |
statistics are calculated independent of each other, memory is |
| 2204 |
– |
allocated for each frame, then freed once correlation calculations are |
| 2205 |
– |
complete for the snapshot. To prevent multiple passes through a |
| 2206 |
– |
potentially large file, \texttt{staticProps} is capable of calculating |
| 2207 |
– |
all requested correlations per frame with only a single pair loop in |
| 2208 |
– |
each frame and a single read of the file. |
| 2209 |
– |
|
| 2210 |
– |
The second, more advanced memory scheme, is used by |
| 2211 |
– |
\texttt{dynamicProps}. Here, the program must have multiple frames in |
| 2212 |
– |
memory to calculate time dependent correlations. In order to prevent a |
| 2213 |
– |
situation where the program runs out of memory due to large |
| 2214 |
– |
trajectories, the user is able to specify that the trajectory be read |
| 2215 |
– |
in blocks. The number of frames in each block is specified by the |
| 2216 |
– |
user, and upon reading a block of the trajectory, |
| 2217 |
– |
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 2218 |
– |
pairs within the block. After in-block correlations are complete, a |
| 2219 |
– |
second block of the trajectory is read, and the cross correlations are |
| 2220 |
– |
calculated between the two blocks. This second block is then freed and |
| 2221 |
– |
then incremented and the process repeated until the end of the |
| 2222 |
– |
trajectory. Once the end is reached, the first block is freed then |
| 2223 |
– |
incremented, and the again the internal time correlations are |
| 2224 |
– |
calculated. The algorithm with the second block is then repeated with |
| 2225 |
– |
the new origin block, until all frame pairs have been correlated in |
| 2226 |
– |
time. This process is illustrated in |
| 2227 |
– |
Fig.~\ref{oopseFig:dynamicPropsMemory}. |
| 2228 |
– |
|
| 2229 |
– |
%\begin{figure} |
| 2230 |
– |
%\centering |
| 2231 |
– |
%\includegraphics[width=\linewidth]{dynamicPropsMem.eps} |
| 2232 |
– |
%\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.} |
| 2233 |
– |
%\label{oopseFig:dynamicPropsMemory} |
| 2234 |
– |
%\end{figure} |
| 2235 |
– |
|
| 1902 |
|
\section{\label{oopseSec:conclusion}Conclusion} |
| 1903 |
|
|
| 1904 |
|
We have presented the design and implementation of our open source |
| 1920 |
|
|
| 1921 |
|
\newpage |
| 1922 |
|
\section{Acknowledgments} |
| 1923 |
< |
The authors would like to thank espresso for fueling this work, and |
| 2258 |
< |
would also like to send a special acknowledgement to single malt |
| 2259 |
< |
scotch for its wonderful calming effects and its ability to make the |
| 2260 |
< |
troubles of the world float away. |
| 1923 |
> |
The authors would like to thank the Notre Dame BoB computer cluster where much of this project was tested. Additionally, the authors would like to acknowledge their funding from {\LARGE FIX ME}. |
| 1924 |
|
|
| 1925 |
|
\bibliographystyle{achemso} |
| 1926 |
|
\bibliography{oopsePaper} |
