| 1 |
|
|
| 2 |
|
\section{\label{sec:mechanics}Mechanics} |
| 3 |
|
|
| 4 |
< |
\subsection{\label{integrate}Integrating the Equations of Motion: the Symplectic Step Integrator} |
| 4 |
> |
\subsection{\label{integrate}Integrating the Equations of Motion: the |
| 5 |
> |
DLM method} |
| 6 |
|
|
| 7 |
< |
Integration of the equations of motion was carried out using the |
| 8 |
< |
symplectic splitting method proposed by Dullweber \emph{et |
| 9 |
< |
al.}.\cite{Dullweber1997} The reason for this integrator selection |
| 10 |
< |
deals with poor energy conservation of rigid body systems using |
| 11 |
< |
quaternions. While quaternions work well for orientational motion in |
| 12 |
< |
alternate ensembles, the microcanonical ensemble has a constant energy |
| 13 |
< |
requirement that is quite sensitive to errors in the equations of |
| 13 |
< |
motion. The original implementation of this code utilized quaternions |
| 14 |
< |
for rotational motion propagation; however, a detailed investigation |
| 15 |
< |
showed that they resulted in a steady drift in the total energy, |
| 16 |
< |
something that has been observed by others.\cite{Laird97} |
| 7 |
> |
The default method for integrating the equations of motion in {\sc |
| 8 |
> |
oopse} is carried out using a velocity-Verlet version of the |
| 9 |
> |
symplectic splitting method proposed by Dullweber, Leimkuhler and |
| 10 |
> |
McLachlan (DLM).\cite{Dullweber1997} When there are no directional |
| 11 |
> |
atoms or rigid bodies present in the simulation, this integrator |
| 12 |
> |
becomes the standard velocity-Verlet integrator which is known to |
| 13 |
> |
sample the microcanonical (NVE) ensemble. |
| 14 |
|
|
| 15 |
+ |
Previous integration methods for orientational motion have problems |
| 16 |
+ |
that are avoided in the DLM method. Direct propagation of the Euler |
| 17 |
+ |
angles has a known $1/\sin\theta$ divergence in the equations of |
| 18 |
+ |
motion for $\phi$ and $\psi$,\cite{AllenTildesley} leading to |
| 19 |
+ |
numerical instabilities any time one of the directional atoms or rigid |
| 20 |
+ |
bodies has an orientation near $\theta=0$ or $\theta=\pi$. More |
| 21 |
+ |
modern quaternion-based integration methods have relatively poor |
| 22 |
+ |
energy conservation. While quaternions work well for orientational |
| 23 |
+ |
motion in other ensembles ensembles, the microcanonical ensemble has a |
| 24 |
+ |
constant energy requirement that is quite sensitive to errors in the |
| 25 |
+ |
equations of motion. An earlier implementation of {\sc oopse} |
| 26 |
+ |
utilized quaternions for propagation of rotational motion; however, a |
| 27 |
+ |
detailed investigation showed that they resulted in a steady drift in |
| 28 |
+ |
the total energy, something that has been observed by |
| 29 |
+ |
others.\cite{Laird97} |
| 30 |
+ |
|
| 31 |
|
The key difference in the integration method proposed by Dullweber |
| 32 |
|
\emph{et al.} is that the entire rotation matrix is propagated from |
| 33 |
< |
one time step to the next. In the past, this would not have been as |
| 34 |
< |
feasible a option, being that the rotation matrix for a single body is |
| 35 |
< |
nine elements long as opposed to 3 or 4 elements for Euler angles and |
| 36 |
< |
quaternions respectively. System memory has become much less of an |
| 37 |
< |
issue in recent times, and this has resulted in substantial benefits |
| 38 |
< |
in energy conservation. There is still the issue of 5 or 6 additional |
| 26 |
< |
elements for describing the orientation of each particle, which will |
| 27 |
< |
increase dump files substantially. Simply translating the rotation |
| 28 |
< |
matrix into its component Euler angles or quaternions for storage |
| 29 |
< |
purposes relieves this burden. |
| 33 |
> |
one time step to the next. In the past, this would not have been a |
| 34 |
> |
feasible option, since that the rotation matrix for a single body has |
| 35 |
> |
nine elements compared to 3 or 4 elements for Euler angles and |
| 36 |
> |
quaternions respectively. System memory has become less costly in |
| 37 |
> |
recent times, and this can be translated into substantial benefits in |
| 38 |
> |
energy conservation. |
| 39 |
|
|
| 40 |
< |
The symplectic splitting method allows for Verlet style integration of |
| 41 |
< |
both linear and angular motion of rigid bodies. In the integration |
| 42 |
< |
method, the orientational propagation involves a sequence of matrix |
| 43 |
< |
evaluations to update the rotation matrix.\cite{Dullweber1997} These |
| 44 |
< |
matrix rotations end up being more costly computationally than the |
| 45 |
< |
simpler arithmetic quaternion propagation. With the same time step, a |
| 46 |
< |
1000 SSD particle simulation shows an average 7\% increase in |
| 47 |
< |
computation time using the symplectic step method in place of |
| 48 |
< |
quaternions. This cost is more than justified when comparing the |
| 49 |
< |
energy conservation of the two methods as illustrated in figure |
| 40 |
> |
The basic equations of motion being integrated are derived from the |
| 41 |
> |
Hamiltonian for conservative systems containing rigid bodies, |
| 42 |
> |
\begin{equation} |
| 43 |
> |
H = \sum_{i} \frac{\vec{v}_i^T M_i \vec{v}_i}{2} + \sum_i |
| 44 |
> |
\frac{\vec{j}_i^T I_i \vec{j_i}}{2} + |
| 45 |
> |
V(\left\{\vec{r}_i\right\}, \left\{\vec{\Omega}_i\right\}) |
| 46 |
> |
\end{equation} |
| 47 |
> |
where $\vec{r}_i$, $\vec{v}_i$, $\vec{j}_i$ and $I_i$ are the |
| 48 |
> |
cartesian position vector of the center of mass, its velocity, the |
| 49 |
> |
body-frame angular velocity, and the moment of inertia tensor (also in |
| 50 |
> |
the body frame) of particle $i$, respectively. $V$ is the potential |
| 51 |
> |
energy function and $\vec{\Omega}_i$ is a representation of the |
| 52 |
> |
orientation of particle $i$. The equations of motion for the particle |
| 53 |
> |
centers of mass are derived from Hamilton's equations and are quite |
| 54 |
> |
simple, |
| 55 |
> |
\begin{eqnarray} |
| 56 |
> |
\frac{d \vec{r}}{d t} & = & \vec{v}(t) \\ |
| 57 |
> |
\frac{d \vec{v}}{d t} & = & \frac{\vec{f}(t)}{m}, |
| 58 |
> |
\end{eqnarray} |
| 59 |
> |
where $\vec{f}(t)$ is the instantaneous force on the centers of mass, |
| 60 |
> |
\begin{equation} |
| 61 |
> |
\vec{f}(t) = \frac{\partial}{\partial |
| 62 |
> |
\vec{r}}V(\left\{\vec{r}_i(t)\right\}, |
| 63 |
> |
\left\{\vec{\Omega}_i(t)\right\}). |
| 64 |
> |
\end{equation} |
| 65 |
> |
|
| 66 |
> |
The equations of motion for the orientational degrees of freedom |
| 67 |
> |
require switching between the space-fixed (or laboratory-fixed) |
| 68 |
> |
frame and the body-fixed frame. Forces and torques are most |
| 69 |
> |
conveniently calculated in the space-fixed frame, while the rotational |
| 70 |
> |
equations of motion are simplest in the body-fixed frame, |
| 71 |
> |
\begin{eqnarray} |
| 72 |
> |
\frac{d j_{x}}{d t} & = & \frac{\tau_x(t)}{I_{xx}} + |
| 73 |
> |
\left(\frac{I_{yy} - I_{zz}}{I_{xx}} \right) j_y j_z \\ |
| 74 |
> |
\frac{d j_{y}}{d t} & = & \frac{\tau_y(t)}{I_{yy}} + |
| 75 |
> |
\left(\frac{I_{zz} - I_{xx}}{I_{yy}} \right) j_z j_x \\ |
| 76 |
> |
\frac{d j_{z}}{d t} & = & \frac{\tau_z(t)}{I_{zz}} + |
| 77 |
> |
\left(\frac{I_{xx} - I_{yy}}{I_{zz}} \right) j_x j_y |
| 78 |
> |
\end{eqnarray} |
| 79 |
> |
|
| 80 |
> |
The DLM method allows for Verlet style integration of both linear and |
| 81 |
> |
angular motion of rigid bodies. |
| 82 |
> |
|
| 83 |
> |
|
| 84 |
> |
|
| 85 |
> |
|
| 86 |
> |
|
| 87 |
> |
In the integration method, the |
| 88 |
> |
orientational propagation involves a sequence of matrix evaluations to |
| 89 |
> |
update the rotation matrix.\cite{Dullweber1997} These matrix rotations |
| 90 |
> |
end up being more costly computationally than the simpler arithmetic |
| 91 |
> |
quaternion propagation. With the same time step, a 1000 SSD particle |
| 92 |
> |
simulation shows an average 7\% increase in computation time using the |
| 93 |
> |
symplectic step method in place of quaternions. This cost is more than |
| 94 |
> |
justified when comparing the energy conservation of the two methods as |
| 95 |
> |
illustrated in figure |
| 96 |
|
\ref{timestep}. |
| 97 |
|
|
| 98 |
|
\begin{figure} |
| 132 |
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 133 |
|
|
| 134 |
|
|
| 135 |
+ |
|
| 136 |
+ |
|
| 137 |
|
{\sc oopse} implements a |
| 138 |
|
|
| 139 |
|
|
