| 1 |
|
| 2 |
\section{\label{sec:mechanics}Mechanics} |
| 3 |
|
| 4 |
\subsection{\label{integrate}Integrating the Equations of Motion: the |
| 5 |
DLM method} |
| 6 |
|
| 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 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 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} |
| 99 |
\epsfxsize=6in |
| 100 |
\epsfbox{timeStep.epsi} |
| 101 |
\caption{Energy conservation using quaternion based integration versus |
| 102 |
the symplectic step method proposed by Dullweber \emph{et al.} with |
| 103 |
increasing time step. For each time step, the dotted line is total |
| 104 |
energy using the symplectic step integrator, and the solid line comes |
| 105 |
from the quaternion integrator. The larger time step plots are shifted |
| 106 |
up from the true energy baseline for clarity.} |
| 107 |
\label{timestep} |
| 108 |
\end{figure} |
| 109 |
|
| 110 |
In figure \ref{timestep}, the resulting energy drift at various time |
| 111 |
steps for both the symplectic step and quaternion integration schemes |
| 112 |
is compared. All of the 1000 SSD particle simulations started with the |
| 113 |
same configuration, and the only difference was the method for |
| 114 |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 115 |
methods for propagating particle rotation conserve energy fairly well, |
| 116 |
with the quaternion method showing a slight energy drift over time in |
| 117 |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 118 |
energy conservation benefits of the symplectic step method are clearly |
| 119 |
demonstrated. Thus, while maintaining the same degree of energy |
| 120 |
conservation, one can take considerably longer time steps, leading to |
| 121 |
an overall reduction in computation time. |
| 122 |
|
| 123 |
Energy drift in these SSD particle simulations was unnoticeable for |
| 124 |
time steps up to three femtoseconds. A slight energy drift on the |
| 125 |
order of 0.012 kcal/mol per nanosecond was observed at a time step of |
| 126 |
four femtoseconds, and as expected, this drift increases dramatically |
| 127 |
with increasing time step. To insure accuracy in the constant energy |
| 128 |
simulations, time steps were set at 2 fs and kept at this value for |
| 129 |
constant pressure simulations as well. |
| 130 |
|
| 131 |
|
| 132 |
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 133 |
|
| 134 |
|
| 135 |
|
| 136 |
|
| 137 |
{\sc oopse} implements a |
| 138 |
|
| 139 |
|
| 140 |
\subsubsection{\label{sec:noseHooverThermo}Nose-Hoover Thermostatting} |
| 141 |
|
| 142 |
To mimic the effects of being in a constant temperature ({\sc nvt}) |
| 143 |
ensemble, {\sc oopse} uses the Nose-Hoover extended system |
| 144 |
approach.\cite{Hoover85} In this method, the equations of motion for |
| 145 |
the particle positions and velocities are |
| 146 |
\begin{eqnarray} |
| 147 |
\dot{{\bf r}} & = & {\bf v} \\ |
| 148 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} |
| 149 |
\label{eq:nosehoovereom} |
| 150 |
\end{eqnarray} |
| 151 |
|
| 152 |
$\chi$ is an ``extra'' variable included in the extended system, and |
| 153 |
it is propagated using the first order equation of motion |
| 154 |
\begin{equation} |
| 155 |
\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right) |
| 156 |
\label{eq:nosehooverext} |
| 157 |
\end{equation} |
| 158 |
where $T_{target}$ is the target temperature for the simulation, and |
| 159 |
$\tau_{T}$ is a time constant for the thermostat. |
| 160 |
|
| 161 |
To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;} |
| 162 |
command would be used in the simulation's {\sc bass} file. There is |
| 163 |
some subtlety in choosing values for $\tau_{T}$, and it is usually set |
| 164 |
to values of a few ps. Within a {\sc bass} file, $\tau_{T}$ could be |
| 165 |
set to 1 ps using the {\tt tauThermostat = 1000; } command. |
| 166 |
|
| 167 |
|
| 168 |
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 169 |
|
| 170 |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 171 |
method was developed to investigate the dynamics of ions inside the ion |
| 172 |
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 173 |
from the deviation of the instaneous force from its mean force. |
| 174 |
|
| 175 |
% |
| 176 |
|
| 177 |
\begin{equation} |
| 178 |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
| 179 |
\end{equation} |
| 180 |
|
| 181 |
|
| 182 |
where% |
| 183 |
\begin{equation} |
| 184 |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle |
| 185 |
\end{equation} |
| 186 |
|
| 187 |
|
| 188 |
If the time-dependent friction decay rapidly, static friction coefficient can |
| 189 |
be approximated by% |
| 190 |
|
| 191 |
\begin{equation} |
| 192 |
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
| 193 |
\end{equation} |
| 194 |
|
| 195 |
|
| 196 |
Hence, diffusion constant can be estimated by |
| 197 |
\begin{equation} |
| 198 |
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 199 |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
| 200 |
\end{equation} |
| 201 |
|
| 202 |
|
| 203 |
\bigskip Z-Constraint method, which fixed the z coordinates of the molecules |
| 204 |
with respect to the center of the mass of the system, was proposed to obtain |
| 205 |
the forces required in force auto-correlation method.\cite{Marrink94} However, |
| 206 |
simply resetting the coordinate will move the center of the mass of the whole |
| 207 |
system. To avoid this problem, a new method was used at {\sc oopse}. Instead of |
| 208 |
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 209 |
well as subtract the total constraint forces from the rest of the system after |
| 210 |
force calculation at each time step. |
| 211 |
\begin{verbatim} |
| 212 |
$F_{\alpha i}=0$ |
| 213 |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$ |
| 214 |
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$ |
| 215 |
$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}}}$ |
| 216 |
\end{verbatim} |
| 217 |
|
| 218 |
At the very beginning of the simulation, the molecules may not be at its |
| 219 |
constraint position. To move the z-constraint molecule to the specified |
| 220 |
position, a simple harmonic potential is used% |
| 221 |
|
| 222 |
\begin{equation} |
| 223 |
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}% |
| 224 |
\end{equation} |
| 225 |
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is |
| 226 |
current z coordinate of the center of mass of the z-constraint molecule, and |
| 227 |
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated |
| 228 |
on the z-constraint molecule at time $t$ can be calculated by% |
| 229 |
\begin{equation} |
| 230 |
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}% |
| 231 |
(z(t)-z_{cons}) |
| 232 |
\end{equation} |
| 233 |
Worthy of mention, other kinds of potential functions can also be used to |
| 234 |
drive the z-constraint molecule. |
