| 1 |
mmeineke |
899 |
|
| 2 |
|
|
\section{\label{sec:mechanics}Mechanics} |
| 3 |
|
|
|
| 4 |
gezelter |
1069 |
\subsection{\label{integrate}Integrating the Equations of Motion: the |
| 5 |
|
|
DLM method} |
| 6 |
mmeineke |
899 |
|
| 7 |
gezelter |
1069 |
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 |
mmeineke |
899 |
|
| 15 |
gezelter |
1069 |
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 |
mmeineke |
899 |
The key difference in the integration method proposed by Dullweber |
| 32 |
|
|
\emph{et al.} is that the entire rotation matrix is propagated from |
| 33 |
gezelter |
1069 |
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 |
mmeineke |
899 |
|
| 40 |
gezelter |
1069 |
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 |
mmeineke |
899 |
\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 |
gezelter |
1069 |
|
| 136 |
|
|
|
| 137 |
mmeineke |
899 |
{\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. |
