| 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 a velocity-Verlet version of the symplectic splitting method |
| 9 |
proposed by Dullweber, Leimkuhler and McLachlan |
| 10 |
(DLM).\cite{Dullweber1997} When there are no directional atoms or |
| 11 |
rigid bodies present in the simulation, this integrator becomes the |
| 12 |
standard velocity-Verlet integrator which is known to sample the |
| 13 |
microcanonical (NVE) ensemble.\cite{} |
| 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, 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 |
| 34 |
feasible, since that the rotation matrix for a single body has nine |
| 35 |
elements compared to 3 or 4 elements for Euler angles and quaternions |
| 36 |
respectively. System memory has become less costly in recent times, |
| 37 |
and this can be translated into substantial benefits in energy |
| 38 |
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} \left( \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 44 |
\frac{1}{2} {\bf j}_i^T \cdot \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot |
| 45 |
{\bf j}_i \right) + |
| 46 |
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right) |
| 47 |
\end{equation} |
| 48 |
where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector |
| 49 |
and velocity of the center of mass of particle $i$, and ${\bf j}_i$ |
| 50 |
and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
| 51 |
momentum and moment of inertia tensor, respectively. $\mathsf{A}_i$ |
| 52 |
is the $3 \times 3$ rotation matrix describing the instantaneous |
| 53 |
orientation of the particle. $V$ is the potential energy function |
| 54 |
which may depend on both the positions $\left\{{\bf r}\right\}$ and |
| 55 |
orientations $\left\{\mathsf{A}\right\}$ of all particles. The |
| 56 |
equations of motion for the particle centers of mass are derived from |
| 57 |
Hamilton's equations and are quite simple, |
| 58 |
\begin{eqnarray} |
| 59 |
\dot{{\bf r}} & = & {\bf v} \\ |
| 60 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} |
| 61 |
\end{eqnarray} |
| 62 |
where ${\bf f}$ is the instantaneous force on the center of mass |
| 63 |
of the particle, |
| 64 |
\begin{equation} |
| 65 |
{\bf f} = - \frac{\partial}{\partial |
| 66 |
{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}). |
| 67 |
\end{equation} |
| 68 |
|
| 69 |
The equations of motion for the orientational degrees of freedom are |
| 70 |
\begin{eqnarray} |
| 71 |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 72 |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 73 |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 74 |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 75 |
V}{\partial \mathsf{A}} \right) |
| 76 |
\end{eqnarray} |
| 77 |
In these equations of motion, the $\mbox{skew}$ matrix of a vector |
| 78 |
${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined: |
| 79 |
\begin{equation} |
| 80 |
\mbox{skew}\left( {\bf v} \right) := \left( |
| 81 |
\begin{array}{ccc} |
| 82 |
0 & v_3 & - v_2 \\ |
| 83 |
-v_3 & 0 & v_1 \\ |
| 84 |
v_2 & -v_1 & 0 |
| 85 |
\end{array} |
| 86 |
\right) |
| 87 |
\end{equation} |
| 88 |
The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$ |
| 89 |
rotation matrix to a vector of orientations by first computing the |
| 90 |
skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and |
| 91 |
then associating this with a length 3 vector by inverting the |
| 92 |
$\mbox{skew}$ function above: |
| 93 |
\begin{equation} |
| 94 |
\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A} |
| 95 |
- \mathsf{A}^{T} \right) |
| 96 |
\end{equation} |
| 97 |
Written this way, the $\mbox{rot}$ operation creates a set of |
| 98 |
conjugate angle coordinates to the body-fixed angular momenta |
| 99 |
represented by ${\bf j}$. This equation of motion for angular momenta |
| 100 |
is equivalent to the more familiar body-fixed forms, |
| 101 |
\begin{eqnarray} |
| 102 |
\dot{j_{x}} & = & \tau^b_x(t) + |
| 103 |
\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z \\ |
| 104 |
\dot{j_{y}} & = & \tau^b_y(t) + |
| 105 |
\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x \\ |
| 106 |
\dot{j_{z}} & = & \tau^b_z(t) + |
| 107 |
\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y |
| 108 |
\end{eqnarray} |
| 109 |
which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are |
| 110 |
most easily derived in the space-fixed frame, |
| 111 |
\begin{equation} |
| 112 |
{\bf \tau}^b(t) = \mathsf{A}(t) \cdot {\bf \tau}^s(t) |
| 113 |
\end{equation} |
| 114 |
where the torques are either derived from the forces on the |
| 115 |
constituent atoms of the rigid body, or for directional atoms, |
| 116 |
directly from derivatives of the potential energy, |
| 117 |
\begin{equation} |
| 118 |
{\bf \tau}^s(t) = - \hat{\bf u}(t) \times \left( \frac{\partial} |
| 119 |
{\partial \hat{\bf u}} V\left(\left\{ {\bf r}(t) \right\}, \left\{ |
| 120 |
\mathsf{A}(t) \right\}\right) \right). |
| 121 |
\end{equation} |
| 122 |
Here $\hat{\bf u}$ is a unit vector pointing along the principal axis |
| 123 |
of the particle in the space-fixed frame. |
| 124 |
|
| 125 |
The DLM method uses a Trotter factorization of the orientational |
| 126 |
propagator. This has three effects: |
| 127 |
\begin{enumerate} |
| 128 |
\item the integrator is area-preserving in phase space (i.e. it is |
| 129 |
{\it symplectic}), |
| 130 |
\item the integrator is time-{\it reversible}, making it suitable for Hybrid |
| 131 |
Monte Carlo applications, and |
| 132 |
\item the error for a single time step is of order $O\left(h^3\right)$ |
| 133 |
for timesteps of length $h$. |
| 134 |
\end{enumerate} |
| 135 |
|
| 136 |
The integration of the equations of motion is carried out in a |
| 137 |
velocity-Verlet style 2 part algorithm: |
| 138 |
|
| 139 |
{\tt moveA:} |
| 140 |
\begin{eqnarray} |
| 141 |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 142 |
v}(t) + \frac{\delta t}{2} \left( {\bf f}(t) / m \right) \\ |
| 143 |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf |
| 144 |
v}\left(t + \delta t / 2 \right) \\ |
| 145 |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 146 |
j}(t) + \frac{\delta t}{2} {\bf \tau}^b(t) \\ |
| 147 |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left( \delta t |
| 148 |
{\bf j}(t + \delta t / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} |
| 149 |
\right) |
| 150 |
\end{eqnarray} |
| 151 |
|
| 152 |
In this context, the $\mathrm{rot}$ function is the reversible product |
| 153 |
of the three body-fixed rotations, |
| 154 |
\begin{equation} |
| 155 |
\mathrm{rot}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
| 156 |
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y / |
| 157 |
2) \cdot \mathsf{G}_x(a_x /2) |
| 158 |
\end{equation} |
| 159 |
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates |
| 160 |
both the rotation matrix ($\mathsf{A}$) and the body-fixed angular |
| 161 |
momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis |
| 162 |
$\alpha$, |
| 163 |
\begin{equation} |
| 164 |
\mathsf{G}_\alpha( \theta ) = \left\{ |
| 165 |
\begin{array}{lcl} |
| 166 |
\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T \\ |
| 167 |
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0) |
| 168 |
\end{array} |
| 169 |
\right. |
| 170 |
\end{equation} |
| 171 |
$\mathsf{R}_\alpha$ is a quadratic approximation to |
| 172 |
the single-axis rotation matrix. For example, in the small-angle |
| 173 |
limit, the rotation matrix around the body-fixed x-axis can be |
| 174 |
approximated as |
| 175 |
\begin{equation} |
| 176 |
\mathsf{R}_x(\theta) = \left( |
| 177 |
\begin{array}{ccc} |
| 178 |
1 & 0 & 0 \\ |
| 179 |
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+ |
| 180 |
\theta^2 / 4} \\ |
| 181 |
0 & \frac{\theta}{1+ |
| 182 |
\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} |
| 183 |
\end{array} |
| 184 |
\right). |
| 185 |
\end{equation} |
| 186 |
All other rotations follow in a straightforward manner. |
| 187 |
|
| 188 |
After the first part of the propagation, the forces and body-fixed |
| 189 |
torques are calculated at the new positions and orientations |
| 190 |
|
| 191 |
{\tt doForces:} |
| 192 |
\begin{eqnarray} |
| 193 |
{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf |
| 194 |
r}(t + \delta t)} \\ |
| 195 |
{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t) |
| 196 |
\times \frac{\partial V}{\partial {\bf u}} \\ |
| 197 |
{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t) |
| 198 |
\cdot {\bf \tau}^s(t + \delta t) |
| 199 |
\end{eqnarray} |
| 200 |
|
| 201 |
{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix |
| 202 |
$\mathsf{A}$ is calculated in {\tt moveA}. Once the forces and |
| 203 |
torques have been obtained at the new time step, the velocities can be |
| 204 |
advanced to the same time value. |
| 205 |
|
| 206 |
{\tt moveB:} |
| 207 |
\begin{eqnarray} |
| 208 |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 209 |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 210 |
{\bf f}(t + \delta t) / m \right) \\ |
| 211 |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 212 |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} {\bf |
| 213 |
\tau}^b(t + \delta t) |
| 214 |
\end{eqnarray} |
| 215 |
|
| 216 |
The matrix rotations used in the DLM method end up being more costly |
| 217 |
computationally than the simpler arithmetic quaternion |
| 218 |
propagation. With the same time step, a 1000-molecule water simulation |
| 219 |
shows an average 7\% increase in computation time using the DLM method |
| 220 |
in place of quaternions. This cost is more than justified when |
| 221 |
comparing the energy conservation of the two methods as illustrated in |
| 222 |
figure \ref{timestep}. |
| 223 |
|
| 224 |
\begin{figure} |
| 225 |
\epsfxsize=6in |
| 226 |
\epsfbox{timeStep.epsi} |
| 227 |
\caption{Energy conservation using quaternion based integration versus |
| 228 |
the method proposed by Dullweber \emph{et al.} with increasing time |
| 229 |
step. For each time step, the dotted line is total energy using the |
| 230 |
symplectic step integrator, and the solid line comes from the |
| 231 |
quaternion integrator. The larger time step plots are shifted up from |
| 232 |
the true energy baseline for clarity.} |
| 233 |
\label{timestep} |
| 234 |
\end{figure} |
| 235 |
|
| 236 |
In figure \ref{timestep}, the resulting energy drift at various time |
| 237 |
steps for both the DLM and quaternion integration schemes is |
| 238 |
compared. All of the 1000 molecule water simulations started with the |
| 239 |
same configuration, and the only difference was the method for |
| 240 |
handling rotational motion. At time steps of 0.1 and 0.5 fs, both |
| 241 |
methods for propagating molecule rotation conserve energy fairly well, |
| 242 |
with the quaternion method showing a slight energy drift over time in |
| 243 |
the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the |
| 244 |
energy conservation benefits of the DLM method are clearly |
| 245 |
demonstrated. Thus, while maintaining the same degree of energy |
| 246 |
conservation, one can take considerably longer time steps, leading to |
| 247 |
an overall reduction in computation time. |
| 248 |
|
| 249 |
|
| 250 |
There is only one specific keyword relevant to the default integrator, |
| 251 |
and that is the time step for integrating the equations of motion. |
| 252 |
|
| 253 |
\begin{tabular}{llll} |
| 254 |
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 255 |
default value} \\ |
| 256 |
$\delta t$ & {\tt dt = 2.0;} & fs & none |
| 257 |
\end{tabular} |
| 258 |
|
| 259 |
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 260 |
|
| 261 |
{\sc oopse} implements a number of extended system integrators for |
| 262 |
sampling from other ensembles relevant to chemical physics. The |
| 263 |
integrator can selected with the {\tt ensemble} keyword in the |
| 264 |
{\tt .bass} file: |
| 265 |
|
| 266 |
\begin{center} |
| 267 |
\begin{tabular}{lll} |
| 268 |
{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\ |
| 269 |
NVE & microcanonical & {\tt ensemble = ``NVE''; } \\ |
| 270 |
NVT & canonical & {\tt ensemble = ``NVT''; } \\ |
| 271 |
NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt |
| 272 |
ensemble = ``NPTi'';} \\ |
| 273 |
NPTf & isobaric-isothermal (with changes to box shape) & {\tt |
| 274 |
ensemble = ``NPTf'';} \\ |
| 275 |
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = |
| 276 |
``NPTxyz'';} \\ |
| 277 |
& (with separate barostats on each box dimension) & |
| 278 |
\end{tabular} |
| 279 |
\end{center} |
| 280 |
|
| 281 |
The relatively well-known Nos\'e-Hoover thermostat is implemented in |
| 282 |
{\sc oopse}'s NVT integrator. This method couples an extra degree of |
| 283 |
freedom (the thermostat) to the kinetic energy of the system, and has |
| 284 |
been shown to sample the canonical distribution in the system degrees |
| 285 |
of freedom while conserving a quantity that is, to within a constant, |
| 286 |
the Helmholtz free energy. |
| 287 |
|
| 288 |
NPT algorithms attempt to maintain constant pressure in the system by |
| 289 |
coupling the volume of the system to a barostat. {\sc oopse} contains |
| 290 |
three different constant pressure algorithms. The first two, NPTi and |
| 291 |
NPTf have been shown to conserve a quantity that is, to within a |
| 292 |
constant, the Gibbs free energy. The Melchionna modification to the |
| 293 |
Hoover barostat is implemented in both NPTi and NPTf. NPTi allows |
| 294 |
only isotropic changes in the simulation box, while box {\it shape} |
| 295 |
variations are allowed in NPTf. The NPTxyz integrator has {\it not} |
| 296 |
been shown to sample from the isobaric-isothermal ensemble. It is |
| 297 |
useful, however, in that it maintains orthogonality for the axes of |
| 298 |
the simulation box while attempting to equalize pressure along the |
| 299 |
three perpendicular directions in the box. |
| 300 |
|
| 301 |
Each of the extended system integrators requires additional keywords |
| 302 |
to set target values for the thermodynamic state variables that are |
| 303 |
being held constant. Keywords are also required to set the |
| 304 |
characteristic decay times for the dynamics of the extended |
| 305 |
variables. |
| 306 |
|
| 307 |
\begin{tabular}{llll} |
| 308 |
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 309 |
default value} \\ |
| 310 |
$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} & K & none \\ |
| 311 |
$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\ |
| 312 |
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\ |
| 313 |
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\ |
| 314 |
& {\tt resetTime = 200;} & fs & none \\ |
| 315 |
& {\tt useInitialExtendedSystemState = ``true'';} & logical & |
| 316 |
false |
| 317 |
\end{tabular} |
| 318 |
|
| 319 |
Two additional keywords can be used to either clear the extended |
| 320 |
system variables periodically ({\tt resetTime}), or to maintain the |
| 321 |
state of the extended system variables between simulations ({\tt |
| 322 |
useInitialExtendedSystemState}). More details on these variables |
| 323 |
and their use in the integrators follows below. |
| 324 |
|
| 325 |
\subsubsection{\label{sec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting} |
| 326 |
|
| 327 |
The Nos\'e-Hoover equations of motion are given by\cite{Hoover85} |
| 328 |
\begin{eqnarray} |
| 329 |
\dot{{\bf r}} & = & {\bf v} \\ |
| 330 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\ |
| 331 |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 332 |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 333 |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 334 |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 335 |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} |
| 336 |
\label{eq:nosehoovereom} |
| 337 |
\end{eqnarray} |
| 338 |
|
| 339 |
$\chi$ is an ``extra'' variable included in the extended system, and |
| 340 |
it is propagated using the first order equation of motion |
| 341 |
\begin{equation} |
| 342 |
\dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right). |
| 343 |
\label{eq:nosehooverext} |
| 344 |
\end{equation} |
| 345 |
|
| 346 |
The instantaneous temperature $T$ is proportional to the total kinetic |
| 347 |
energy (both translational and orientational) and is given by |
| 348 |
\begin{equation} |
| 349 |
T = \frac{2 K}{f k_B} |
| 350 |
\end{equation} |
| 351 |
Here, $f$ is the total number of degrees of freedom in the system, |
| 352 |
\begin{equation} |
| 353 |
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}} |
| 354 |
\end{equation} |
| 355 |
and $K$ is the total kinetic energy, |
| 356 |
\begin{equation} |
| 357 |
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 358 |
\sum_{i=1}^{N_{\mathrm{orient}}} \sum_{\alpha=x,y,z} \frac{{\bf |
| 359 |
j}_{i\alpha}^T \cdot {\bf j}_{i\alpha}}{2 |
| 360 |
\overleftrightarrow{\mathsf{I}}_{i,\alpha \alpha}} |
| 361 |
\end{equation} |
| 362 |
|
| 363 |
In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for |
| 364 |
relaxation of the temperature to the target value. To set values for |
| 365 |
$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the |
| 366 |
{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt |
| 367 |
.bass} file. The units for {\tt tauThermostat} are fs, and the units |
| 368 |
for the {\tt targetTemperature} are degrees K. The integration of |
| 369 |
the equations of motion is carried out in a velocity-Verlet style 2 |
| 370 |
part algorithm: |
| 371 |
|
| 372 |
{\tt moveA:} |
| 373 |
\begin{eqnarray} |
| 374 |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 375 |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 376 |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 377 |
\chi(t)\right) \\ |
| 378 |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf |
| 379 |
v}\left(t + \delta t / 2 \right) \\ |
| 380 |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 381 |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 382 |
\chi(t) \right) \\ |
| 383 |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 384 |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{b} \right) \\ |
| 385 |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 386 |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 387 |
\right) |
| 388 |
\end{eqnarray} |
| 389 |
|
| 390 |
Here $\mathrm{rot}( {\bf j} / {\bf I} )$ is the same symplectic |
| 391 |
Trotter factorization of the three rotation operations that was |
| 392 |
discussed in the section on the DLM integrator. Note that this |
| 393 |
operation modifies both the rotation matrix $\mathsf{A}$ and the |
| 394 |
angular momentum ${\bf j}$. {\tt moveA} propagates velocities by a |
| 395 |
half time step, and positional degrees of freedom by a full time step. |
| 396 |
The new positions (and orientations) are then used to calculate a new |
| 397 |
set of forces and torques. |
| 398 |
|
| 399 |
{\tt doForces:} |
| 400 |
\begin{eqnarray} |
| 401 |
{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf |
| 402 |
r}(t + \delta t)} \\ |
| 403 |
{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t) |
| 404 |
\times \frac{\partial V}{\partial {\bf u}} \\ |
| 405 |
{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t) |
| 406 |
\cdot {\bf \tau}^s(t + \delta t) |
| 407 |
\end{eqnarray} |
| 408 |
|
| 409 |
Here ${\bf u}$ is a unit vector pointing along the principal axis of |
| 410 |
the rigid body being propagated, ${\bf \tau}^s$ is the torque in the |
| 411 |
space-fixed (laboratory) frame, and ${\bf \tau}^b$ is the torque in |
| 412 |
the body-fixed frame. ${\bf u}$ is automatically calculated when the |
| 413 |
rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}. |
| 414 |
|
| 415 |
Once the forces and torques have been obtained at the new time step, |
| 416 |
the velocities can be advanced to the same time value. |
| 417 |
|
| 418 |
{\tt moveB:} |
| 419 |
\begin{eqnarray} |
| 420 |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 421 |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 422 |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 423 |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 424 |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 425 |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 426 |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 427 |
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t) |
| 428 |
\chi(t \delta t)\right) \\ |
| 429 |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 430 |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 431 |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 432 |
\chi(t + \delta t) \right) |
| 433 |
\end{eqnarray} |
| 434 |
|
| 435 |
Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required |
| 436 |
to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, they |
| 437 |
indirectly depend on their own values at time $t + \delta t$. {\tt |
| 438 |
moveB} is therefore done in an iterative fashion until $\chi(t + |
| 439 |
\delta t)$ becomes self-consistent. The relative tolerance for the |
| 440 |
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but |
| 441 |
{\sc oopse} will terminate the iteration after 4 loops even if the |
| 442 |
consistency check has not been satisfied. |
| 443 |
|
| 444 |
The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the |
| 445 |
extended system that is, to within a constant, identical to the |
| 446 |
Helmholtz free energy, |
| 447 |
\begin{equation} |
| 448 |
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 449 |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime |
| 450 |
\right) |
| 451 |
\end{equation} |
| 452 |
Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation |
| 453 |
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the |
| 454 |
last column of the {\tt .stat} file to allow checks on the quality of |
| 455 |
the integration. |
| 456 |
|
| 457 |
Bond constraints are applied at the end of both the {\tt moveA} and |
| 458 |
{\tt moveB} portions of the algorithm. Details on the constraint |
| 459 |
algorithms are given in section \ref{sec:rattle}. |
| 460 |
|
| 461 |
\subsubsection{\label{sec:NPTi}Constant-pressure integration (isotropic box)} |
| 462 |
|
| 463 |
To carry out isobaric-isothermal ensemble calculations {\sc oopse} |
| 464 |
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen |
| 465 |
equations of motion,\cite{Melchionna93} |
| 466 |
|
| 467 |
\begin{eqnarray} |
| 468 |
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right) \\ |
| 469 |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v} \\ |
| 470 |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 471 |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 472 |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 473 |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 474 |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\ |
| 475 |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 476 |
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\ |
| 477 |
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - |
| 478 |
P_{\mathrm{target}} \right) \\ |
| 479 |
\dot{V} & = & 3 V \eta |
| 480 |
\label{eq:melchionna1} |
| 481 |
\end{eqnarray} |
| 482 |
|
| 483 |
$\chi$ and $\eta$ are the ``extra'' degrees of freedom in the extended |
| 484 |
system. $\chi$ is a thermostat, and it has the same function as it |
| 485 |
does in the Nos\'e-Hoover NVT integrator. $\eta$ is a barostat which |
| 486 |
controls changes to the volume of the simulation box. ${\bf R}_0$ is |
| 487 |
the location of the center of mass for the entire system. |
| 488 |
|
| 489 |
The NPTi integrator requires an instantaneous pressure. This quantity |
| 490 |
is calculated via the pressure tensor, |
| 491 |
\begin{equation} |
| 492 |
\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{V(t)} \left( \sum_{i=1}^{N} |
| 493 |
m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) + |
| 494 |
\overleftrightarrow{\mathsf{W}}(t) |
| 495 |
\end{equation} |
| 496 |
The kinetic contribution to the pressure tensor utilizes the {\it |
| 497 |
outer} product of the velocities denoted by the $\otimes$ symbol. The |
| 498 |
stress tensor is calculated from another outer product of the |
| 499 |
inter-atomic separation vectors (${\bf r}_{ij} = {\bf r}_j - {\bf |
| 500 |
r}_i$) with the forces between the same two atoms, |
| 501 |
\begin{equation} |
| 502 |
\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t) |
| 503 |
\otimes {\bf f}_{ij}(t) |
| 504 |
\end{equation} |
| 505 |
The instantaneous pressure is then simply obtained from the trace of |
| 506 |
the Pressure tensor, |
| 507 |
\begin{equation} |
| 508 |
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 509 |
\right) |
| 510 |
\end{equation} |
| 511 |
|
| 512 |
In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for |
| 513 |
relaxation of the pressure to the target value. To set values for |
| 514 |
$\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the |
| 515 |
{\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt .bass} |
| 516 |
file. The units for {\tt tauBarostat} are fs, and the units for the |
| 517 |
{\tt targetPressure} are atmospheres. Like in the NVT integrator, the |
| 518 |
integration of the equations of motion is carried out in a |
| 519 |
velocity-Verlet style 2 part algorithm: |
| 520 |
|
| 521 |
{\tt moveA:} |
| 522 |
\begin{eqnarray} |
| 523 |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 524 |
P(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\ |
| 525 |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 526 |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 527 |
\left(\chi(t) + \eta(t) \right) \right) \\ |
| 528 |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 529 |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 530 |
\chi(t) \right) \\ |
| 531 |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 532 |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{b} \right) \\ |
| 533 |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 534 |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 535 |
\right) \\ |
| 536 |
\eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t V(t)}{2 N k_B |
| 537 |
T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\ |
| 538 |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf |
| 539 |
v}\left(t + \delta t / 2 \right) + \eta(t + \delta t / 2)\left[ {\bf |
| 540 |
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\ |
| 541 |
\mathsf{H}(t + \delta t) & \leftarrow & e^{-\delta t \eta(t + \delta t |
| 542 |
/ 2)} \mathsf{H}(t) |
| 543 |
\end{eqnarray} |
| 544 |
|
| 545 |
Most of these equations are identical to their counterparts in the NVT |
| 546 |
integrator, but the propagation of positions to time $t + \delta t$ |
| 547 |
depends on the positions at the same time. {\sc oopse} carries out |
| 548 |
this step iteratively (with a limit of 5 passes through the iterative |
| 549 |
loop). Also, the simulation box $\mathsf{H}$ is scaled uniformly for |
| 550 |
one full time step by an exponential factor that depends on the value |
| 551 |
of $\eta$ at time $t + |
| 552 |
\delta t / 2$. Reshaping the box uniformly also scales the volume of |
| 553 |
the box by |
| 554 |
\begin{equation} |
| 555 |
V(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)} |
| 556 |
V(t) |
| 557 |
\end{equation} |
| 558 |
|
| 559 |
The {\tt doForces} step for the NPTi integrator is exactly the same as |
| 560 |
in the DLM and NVT integrators. Once the forces and torques have been |
| 561 |
obtained at the new time step, the velocities can be advanced to the |
| 562 |
same time value. |
| 563 |
|
| 564 |
{\tt moveB:} |
| 565 |
\begin{eqnarray} |
| 566 |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 567 |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 568 |
P(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\}, |
| 569 |
\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\ |
| 570 |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 571 |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 572 |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 573 |
\eta(t + \delta t) & \leftarrow & \eta(t + \delta t / 2) + |
| 574 |
\frac{\delta t V(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2} |
| 575 |
\left( P(t + \delta t) - P_{\mathrm{target}} |
| 576 |
\right) \\ |
| 577 |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 578 |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 579 |
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t) |
| 580 |
(\chi(t + \delta t) + \eta(t + \delta t)) \right) \\ |
| 581 |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 582 |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 583 |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 584 |
\chi(t + \delta t) \right) |
| 585 |
\end{eqnarray} |
| 586 |
|
| 587 |
Once again, since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ |
| 588 |
are required to caclculate $T(t + \delta t)$, $P(t + \delta t)$, $\chi(t + |
| 589 |
\delta t)$, and $\eta(t + \delta t)$, they indirectly depend on their |
| 590 |
own values at time $t + \delta t$. {\tt moveB} is therefore done in |
| 591 |
an iterative fashion until $\chi(t + \delta t)$ and $\eta(t + \delta |
| 592 |
t)$ become self-consistent. The relative tolerance for the |
| 593 |
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but |
| 594 |
{\sc oopse} will terminate the iteration after 4 loops even if the |
| 595 |
consistency check has not been satisfied. |
| 596 |
|
| 597 |
The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is |
| 598 |
known to conserve a Hamiltonian for the extended system that is, to |
| 599 |
within a constant, identical to the Gibbs free energy, |
| 600 |
\begin{equation} |
| 601 |
H_{\mathrm{NPTi}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 602 |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime |
| 603 |
\right) + P_{\mathrm{target}} V(t). |
| 604 |
\end{equation} |
| 605 |
Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can result in |
| 606 |
non-conservation of $H_{\mathrm{NPTi}}$, so the conserved quantity is |
| 607 |
maintained in the last column of the {\tt .stat} file to allow checks |
| 608 |
on the quality of the integration. It is also known that this |
| 609 |
algorithm samples the equilibrium distribution for the enthalpy |
| 610 |
(including contributions for the thermostat and barostat), |
| 611 |
\begin{equation} |
| 612 |
H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2} \left( |
| 613 |
\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) + P_{\mathrm{target}} V(t). |
| 614 |
\end{equation} |
| 615 |
|
| 616 |
Bond constraints are applied at the end of both the {\tt moveA} and |
| 617 |
{\tt moveB} portions of the algorithm. Details on the constraint |
| 618 |
algorithms are given in section \ref{sec:rattle}. |
| 619 |
+ |
| 620 |
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 621 |
|
| 622 |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 623 |
method was developed to investigate the dynamics of ions inside the ion |
| 624 |
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 625 |
from the deviation of the instaneous force from its mean force. |
| 626 |
|
| 627 |
% |
| 628 |
|
| 629 |
\begin{equation} |
| 630 |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
| 631 |
\end{equation} |
| 632 |
|
| 633 |
|
| 634 |
where% |
| 635 |
\begin{equation} |
| 636 |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle |
| 637 |
\end{equation} |
| 638 |
|
| 639 |
|
| 640 |
If the time-dependent friction decay rapidly, static friction coefficient can |
| 641 |
be approximated by% |
| 642 |
|
| 643 |
\begin{equation} |
| 644 |
\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
| 645 |
\end{equation} |
| 646 |
|
| 647 |
|
| 648 |
Hence, diffusion constant can be estimated by |
| 649 |
\begin{equation} |
| 650 |
D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 651 |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
| 652 |
\end{equation} |
| 653 |
|
| 654 |
|
| 655 |
\bigskip Z-Constraint method, which fixed the z coordinates of the molecules |
| 656 |
with respect to the center of the mass of the system, was proposed to obtain |
| 657 |
the forces required in force auto-correlation method.\cite{Marrink94} However, |
| 658 |
simply resetting the coordinate will move the center of the mass of the whole |
| 659 |
system. To avoid this problem, a new method was used at {\sc oopse}. Instead of |
| 660 |
resetting the coordinate, we reset the forces of z-constraint molecules as |
| 661 |
well as subtract the total constraint forces from the rest of the system after |
| 662 |
force calculation at each time step. |
| 663 |
\begin{verbatim} |
| 664 |
$F_{\alpha i}=0$ |
| 665 |
$V_{\alpha i}=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}$ |
| 666 |
$F_{\alpha i}=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}$ |
| 667 |
$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}}}$ |
| 668 |
\end{verbatim} |
| 669 |
|
| 670 |
At the very beginning of the simulation, the molecules may not be at its |
| 671 |
constraint position. To move the z-constraint molecule to the specified |
| 672 |
position, a simple harmonic potential is used% |
| 673 |
|
| 674 |
\begin{equation} |
| 675 |
U(t)=\frac{1}{2}k_{Harmonic}(z(t)-z_{cons})^{2}% |
| 676 |
\end{equation} |
| 677 |
where $k_{Harmonic}$\bigskip\ is the harmonic force constant, $z(t)$ is |
| 678 |
current z coordinate of the center of mass of the z-constraint molecule, and |
| 679 |
$z_{cons}$ is the restraint position. Therefore, the harmonic force operated |
| 680 |
on the z-constraint molecule at time $t$ can be calculated by% |
| 681 |
\begin{equation} |
| 682 |
F_{z_{Harmonic}}(t)=-\frac{\partial U(t)}{\partial z(t)}=-k_{Harmonic}% |
| 683 |
(z(t)-z_{cons}) |
| 684 |
\end{equation} |
| 685 |
Worthy of mention, other kinds of potential functions can also be used to |
| 686 |
drive the z-constraint molecule. |
