| 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} |
| 29 |
> |
Laird {\it et al.}\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. |
| 32 |
> |
\emph{et al.} is that the entire $3 \times 3$ rotation matrix is |
| 33 |
> |
propagated from one time step to the next. In the past, this would not |
| 34 |
> |
have been feasible, since the rotation matrix for a single body has |
| 35 |
> |
nine elements compared with the more memory-efficient methods (using |
| 36 |
> |
three Euler angles or 4 quaternions). Computer memory has become much |
| 37 |
> |
less costly in recent years, and this can be translated into |
| 38 |
> |
substantial benefits in energy conservation. |
| 39 |
|
|
| 40 |
|
The basic equations of motion being integrated are derived from the |
| 41 |
|
Hamiltonian for conservative systems containing rigid bodies, |
| 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} |
| 72 |
> |
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\ |
| 73 |
> |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{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} |
| 134 |
|
\end{enumerate} |
| 135 |
|
|
| 136 |
|
The integration of the equations of motion is carried out in a |
| 137 |
< |
velocity-Verlet style 2 part algorithm: |
| 137 |
> |
velocity-Verlet style 2-part algorithm: |
| 138 |
|
|
| 139 |
|
{\tt moveA:} |
| 140 |
|
\begin{eqnarray} |
| 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( |
| 176 |
> |
\mathsf{R}_x(\theta) \approx \left( |
| 177 |
|
\begin{array}{ccc} |
| 178 |
|
1 & 0 & 0 \\ |
| 179 |
|
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+ |
| 190 |
|
|
| 191 |
|
{\tt doForces:} |
| 192 |
|
\begin{eqnarray} |
| 193 |
< |
{\bf f}(t + \delta t) & \leftarrow & - \frac{\partial V}{\partial {\bf |
| 194 |
< |
r}(t + \delta t)} \\ |
| 193 |
> |
{\bf f}(t + \delta t) & \leftarrow & - \left(\frac{\partial V}{\partial {\bf |
| 194 |
> |
r}}\right)_{{\bf 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) |
| 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.} |
| 230 |
> |
DLM integrator, and the solid line comes from the quaternion |
| 231 |
> |
integrator. The larger time step plots are shifted up from the true |
| 232 |
> |
energy baseline for clarity.} |
| 233 |
|
\label{timestep} |
| 234 |
|
\end{figure} |
| 235 |
|
|
| 246 |
|
conservation, one can take considerably longer time steps, leading to |
| 247 |
|
an overall reduction in computation time. |
| 248 |
|
|
| 249 |
– |
|
| 249 |
|
There is only one specific keyword relevant to the default integrator, |
| 250 |
|
and that is the time step for integrating the equations of motion. |
| 251 |
|
|
| 252 |
+ |
\begin{center} |
| 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 |
+ |
\end{center} |
| 259 |
|
|
| 260 |
|
\subsection{\label{sec:extended}Extended Systems for other Ensembles} |
| 261 |
|
|
| 330 |
|
\dot{{\bf r}} & = & {\bf v} \\ |
| 331 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\ |
| 332 |
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 333 |
< |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 334 |
< |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 333 |
> |
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\ |
| 334 |
> |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} |
| 335 |
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 336 |
|
V}{\partial \mathsf{A}} \right) - \chi {\bf j} |
| 337 |
|
\label{eq:nosehoovereom} |
| 356 |
|
and $K$ is the total kinetic energy, |
| 357 |
|
\begin{equation} |
| 358 |
|
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 359 |
< |
\sum_{i=1}^{N_{\mathrm{orient}}} \sum_{\alpha=x,y,z} \frac{{\bf |
| 360 |
< |
j}_{i\alpha}^T \cdot {\bf j}_{i\alpha}}{2 |
| 360 |
< |
\overleftrightarrow{\mathsf{I}}_{i,\alpha \alpha}} |
| 359 |
> |
\sum_{i=1}^{N_{\mathrm{orient}}} \frac{1}{2} {\bf j}_i^T \cdot |
| 360 |
> |
\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i |
| 361 |
|
\end{equation} |
| 362 |
|
|
| 363 |
|
In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for |
| 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) \\ |
| 384 |
> |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \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) |
| 388 |
|
\end{eqnarray} |
| 389 |
|
|
| 390 |
< |
Here ${\bf u}$ is a unit vector pointing along the principal axis of |
| 391 |
< |
the rigid body being propagated, ${\bf \tau}^s$ is the torque in the |
| 392 |
< |
space-fixed (laboratory) frame, and ${\bf \tau}^b$ is the torque in |
| 393 |
< |
the body-fixed frame. ${\bf u}$ is automatically calculated when the |
| 394 |
< |
rotation matrix $\mathsf{A}$ is calculated in {\tt moveA}. |
| 390 |
> |
Here $\mathrm{rot}(\delta t * {\bf j} |
| 391 |
> |
\overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Trotter |
| 392 |
> |
factorization of the three rotation operations that was discussed in |
| 393 |
> |
the section on the DLM integrator. Note that this operation modifies |
| 394 |
> |
both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf |
| 395 |
> |
j}$. {\tt moveA} propagates velocities by a half time step, and |
| 396 |
> |
positional degrees of freedom by a full time step. The new positions |
| 397 |
> |
(and orientations) are then used to calculate a new set of forces and |
| 398 |
> |
torques in exactly the same way they are calculated in the {\tt |
| 399 |
> |
doForces} portion of the DLM integrator. |
| 400 |
|
|
| 401 |
|
Once the forces and torques have been obtained at the new time step, |
| 402 |
< |
the velocities can be advanced to the same time value. |
| 402 |
> |
the temperature, velocities, and the extended system variable can be |
| 403 |
> |
advanced to the same time value. |
| 404 |
|
|
| 405 |
|
{\tt moveB:} |
| 406 |
|
\begin{eqnarray} |
| 433 |
|
Helmholtz free energy, |
| 434 |
|
\begin{equation} |
| 435 |
|
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 436 |
< |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime |
| 436 |
> |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
| 437 |
|
\right) |
| 438 |
|
\end{equation} |
| 439 |
|
Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation |
| 445 |
|
{\tt moveB} portions of the algorithm. Details on the constraint |
| 446 |
|
algorithms are given in section \ref{sec:rattle}. |
| 447 |
|
|
| 448 |
< |
\subsubsection{\label{sec:NPTi}Constant-pressure integration (isotropic box)} |
| 448 |
> |
\subsubsection{\label{sec:NPTi}Constant-pressure integration with |
| 449 |
> |
isotropic box deformations (NPTi)} |
| 450 |
|
|
| 451 |
|
To carry out isobaric-isothermal ensemble calculations {\sc oopse} |
| 452 |
|
implements the Melchionna modifications to the Nos\'e-Hoover-Andersen |
| 464 |
|
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\ |
| 465 |
|
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - |
| 466 |
|
P_{\mathrm{target}} \right) \\ |
| 467 |
< |
\dot{V} & = & 3 V \eta |
| 467 |
> |
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta |
| 468 |
|
\label{eq:melchionna1} |
| 469 |
|
\end{eqnarray} |
| 470 |
|
|
| 472 |
|
system. $\chi$ is a thermostat, and it has the same function as it |
| 473 |
|
does in the Nos\'e-Hoover NVT integrator. $\eta$ is a barostat which |
| 474 |
|
controls changes to the volume of the simulation box. ${\bf R}_0$ is |
| 475 |
< |
the location of the center of mass for the entire system. |
| 475 |
> |
the location of the center of mass for the entire system, and |
| 476 |
> |
$\mathcal{V}$ is the volume of the simulation box. At any time, the |
| 477 |
> |
volume can be calculated from the determinant of the matrix which |
| 478 |
> |
describes the box shape: |
| 479 |
> |
\begin{equation} |
| 480 |
> |
\mathcal{V} = \det(\mathsf{H}) |
| 481 |
> |
\end{equation} |
| 482 |
|
|
| 483 |
|
The NPTi integrator requires an instantaneous pressure. This quantity |
| 484 |
|
is calculated via the pressure tensor, |
| 485 |
|
\begin{equation} |
| 486 |
< |
\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{V(t)} \left( \sum_{i=1}^{N} |
| 487 |
< |
m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) + |
| 486 |
> |
\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{\mathcal{V}(t)} \left( |
| 487 |
> |
\sum_{i=1}^{N} m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) + |
| 488 |
|
\overleftrightarrow{\mathsf{W}}(t) |
| 489 |
|
\end{equation} |
| 490 |
|
The kinetic contribution to the pressure tensor utilizes the {\it |
| 523 |
|
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 524 |
|
\chi(t) \right) \\ |
| 525 |
|
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 526 |
< |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{b} \right) \\ |
| 526 |
> |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\ |
| 527 |
|
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 528 |
|
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 529 |
|
\right) \\ |
| 530 |
< |
\eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t V(t)}{2 N k_B |
| 531 |
< |
T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\ |
| 530 |
> |
\eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B |
| 531 |
> |
T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\ |
| 532 |
|
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf |
| 533 |
|
v}\left(t + \delta t / 2 \right) + \eta(t + \delta t / 2)\left[ {\bf |
| 534 |
|
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\ |
| 546 |
|
\delta t / 2$. Reshaping the box uniformly also scales the volume of |
| 547 |
|
the box by |
| 548 |
|
\begin{equation} |
| 549 |
< |
V(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)} |
| 550 |
< |
V(t) |
| 549 |
> |
\mathcal{V}(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)} |
| 550 |
> |
\mathcal{V}(t) |
| 551 |
|
\end{equation} |
| 552 |
|
|
| 553 |
|
The {\tt doForces} step for the NPTi integrator is exactly the same as |
| 554 |
< |
in the DLM and NVT integrators. Once the forces and torques have been |
| 555 |
< |
obtained at the new time step, the velocities can be advanced to the |
| 556 |
< |
same time value. |
| 554 |
> |
in both the DLM and NVT integrators. Once the forces and torques have |
| 555 |
> |
been obtained at the new time step, the velocities can be advanced to |
| 556 |
> |
the same time value. |
| 557 |
|
|
| 558 |
|
{\tt moveB:} |
| 559 |
|
\begin{eqnarray} |
| 565 |
|
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 566 |
|
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 567 |
|
\eta(t + \delta t) & \leftarrow & \eta(t + \delta t / 2) + |
| 568 |
< |
\frac{\delta t V(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2} |
| 568 |
> |
\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2} |
| 569 |
|
\left( P(t + \delta t) - P_{\mathrm{target}} |
| 570 |
|
\right) \\ |
| 571 |
|
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 593 |
|
within a constant, identical to the Gibbs free energy, |
| 594 |
|
\begin{equation} |
| 595 |
|
H_{\mathrm{NPTi}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 596 |
< |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t\prime) dt\prime |
| 597 |
< |
\right) + P_{\mathrm{target}} V(t). |
| 596 |
> |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
| 597 |
> |
\right) + P_{\mathrm{target}} \mathcal{V}(t). |
| 598 |
|
\end{equation} |
| 599 |
|
Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can result in |
| 600 |
|
non-conservation of $H_{\mathrm{NPTi}}$, so the conserved quantity is |
| 604 |
|
(including contributions for the thermostat and barostat), |
| 605 |
|
\begin{equation} |
| 606 |
|
H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2} \left( |
| 607 |
< |
\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) + P_{\mathrm{target}} V(t). |
| 607 |
> |
\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) + P_{\mathrm{target}} |
| 608 |
> |
\mathcal{V}(t). |
| 609 |
|
\end{equation} |
| 610 |
|
|
| 611 |
|
Bond constraints are applied at the end of both the {\tt moveA} and |
| 612 |
|
{\tt moveB} portions of the algorithm. Details on the constraint |
| 613 |
|
algorithms are given in section \ref{sec:rattle}. |
| 614 |
< |
+ |
| 614 |
> |
|
| 615 |
> |
\subsubsection{\label{sec:NPTf}Constant-pressure integration with a |
| 616 |
> |
flexible box (NPTf)} |
| 617 |
> |
|
| 618 |
> |
There is a relatively simple generalization of the |
| 619 |
> |
Nos\'e-Hoover-Andersen method to include changes in the simulation box |
| 620 |
> |
{\it shape} as well as in the volume of the box. This method utilizes |
| 621 |
> |
the full $3 \times 3$ pressure tensor and introduces a tensor of |
| 622 |
> |
extended variables ($\overleftrightarrow{\eta}$) to control changes to |
| 623 |
> |
the box shape. The equations of motion for this method are |
| 624 |
> |
\begin{eqnarray} |
| 625 |
> |
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right) \\ |
| 626 |
> |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} + |
| 627 |
> |
\chi \mathsf{1}) {\bf v} \\ |
| 628 |
> |
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 629 |
> |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 630 |
> |
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 631 |
> |
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 632 |
> |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\ |
| 633 |
> |
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 634 |
> |
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\ |
| 635 |
> |
\dot{\overleftrightarrow{eta}} & = & \frac{1}{\tau_{B}^2 f k_B |
| 636 |
> |
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) \\ |
| 637 |
> |
\dot{\mathsf{H}} & = & \overleftrightarrow{\eta} \cdot \mathsf{H} |
| 638 |
> |
\label{eq:melchionna2} |
| 639 |
> |
\end{eqnarray} |
| 640 |
> |
|
| 641 |
> |
Here, $\mathsf{1}$ is the unit matrix and $\overleftrightarrow{\mathsf{P}}$ |
| 642 |
> |
is the pressure tensor. Again, the volume, $\mathcal{V} = \det |
| 643 |
> |
\mathsf{H}$. |
| 644 |
> |
|
| 645 |
> |
The propagation of the equations of motion is nearly identical to the |
| 646 |
> |
NPTi integration: |
| 647 |
> |
|
| 648 |
> |
{\tt moveA:} |
| 649 |
> |
\begin{eqnarray} |
| 650 |
> |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 651 |
> |
\overleftrightarrow{\mathsf{P}}(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\ |
| 652 |
> |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 653 |
> |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - |
| 654 |
> |
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
| 655 |
> |
{\bf v}(t) \right) \\ |
| 656 |
> |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 657 |
> |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 658 |
> |
\chi(t) \right) \\ |
| 659 |
> |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 660 |
> |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\ |
| 661 |
> |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 662 |
> |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 663 |
> |
\right) \\ |
| 664 |
> |
\overleftrightarrow{\eta}(t + \delta t / 2) & \leftarrow & \overleftrightarrow{\eta}(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B |
| 665 |
> |
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) - P_{\mathrm{target}}\mathsf{1} \right) \\ |
| 666 |
> |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf |
| 667 |
> |
v}\left(t + \delta t / 2 \right) + \overleftrightarrow{\eta}(t + |
| 668 |
> |
\delta t / 2) \cdot \left[ {\bf |
| 669 |
> |
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\ |
| 670 |
> |
\mathsf{H}(t + \delta t) & \leftarrow & \mathsf{H}(t) \cdot e^{-\delta t |
| 671 |
> |
\overleftrightarrow{\eta}(t + \delta t / 2)} |
| 672 |
> |
\end{eqnarray} |
| 673 |
> |
{\sc oopse} uses a power series expansion truncated at second order |
| 674 |
> |
for the exponential operation which scales the simulation box. |
| 675 |
> |
|
| 676 |
> |
The {\tt moveB} portion of the algorithm is largely unchanged from the |
| 677 |
> |
NPTi integrator: |
| 678 |
> |
|
| 679 |
> |
{\tt moveB:} |
| 680 |
> |
\begin{eqnarray} |
| 681 |
> |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 682 |
> |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 683 |
> |
\overleftrightarrow{\mathsf{P}}(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\}, |
| 684 |
> |
\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\ |
| 685 |
> |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 686 |
> |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 687 |
> |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 688 |
> |
\overleftrightarrow{\eta}(t + \delta t) & \leftarrow & \overleftrightarrow{\eta}(t + \delta t / 2) + |
| 689 |
> |
\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2} |
| 690 |
> |
\left( \overleftrightarrow{P}(t + \delta t) - P_{\mathrm{target}}\mathsf{1} |
| 691 |
> |
\right) \\ |
| 692 |
> |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 693 |
> |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 694 |
> |
\frac{{\bf f}(t + \delta t)}{m} - |
| 695 |
> |
(\chi(t + \delta t)\mathsf{1} + \overleftrightarrow{\eta}(t + \delta |
| 696 |
> |
t)) \right) \cdot {\bf v}(t + \delta t) \\ |
| 697 |
> |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 698 |
> |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 699 |
> |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 700 |
> |
\chi(t + \delta t) \right) |
| 701 |
> |
\end{eqnarray} |
| 702 |
> |
|
| 703 |
> |
The iterative schemes for both {\tt moveA} and {\tt moveB} are |
| 704 |
> |
identical to those described for the NPTi integrator. |
| 705 |
> |
|
| 706 |
> |
The NPTf integrator is known to conserve the following Hamiltonian: |
| 707 |
> |
\begin{equation} |
| 708 |
> |
H_{\mathrm{NPTf}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 709 |
> |
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
| 710 |
> |
\right) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B |
| 711 |
> |
T_{\mathrm{target}}}{2} |
| 712 |
> |
\mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2. |
| 713 |
> |
\end{equation} |
| 714 |
> |
|
| 715 |
> |
This integrator must be used with care, particularly in liquid |
| 716 |
> |
simulations. Liquids have very small restoring forces in the |
| 717 |
> |
off-diagonal directions, and the simulation box can very quickly form |
| 718 |
> |
elongated and sheared geometries which become smaller than the |
| 719 |
> |
electrostatic or Lennard-Jones cutoff radii. It finds most use in |
| 720 |
> |
simulating crystals or liquid crystals which assume non-orthorhombic |
| 721 |
> |
geometries. |
| 722 |
> |
|
| 723 |
> |
\subsubsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)} |
| 724 |
> |
|
| 725 |
> |
There is one additional extended system integrator which is somewhat |
| 726 |
> |
simpler than the NPTf method described above. In this case, the three |
| 727 |
> |
axes have independent barostats which each attempt to preserve the |
| 728 |
> |
target pressure along the box walls perpendicular to that particular |
| 729 |
> |
axis. The lengths of the box axes are allowed to fluctuate |
| 730 |
> |
independently, but the angle between the box axes does not change. |
| 731 |
> |
The equations of motion are identical to those described above, but |
| 732 |
> |
only the {\it diagonal} elements of $\overleftrightarrow{\eta}$ are |
| 733 |
> |
computed. The off-diagonal elements are set to zero (even when the |
| 734 |
> |
pressure tensor has non-zero off-diagonal elements). |
| 735 |
> |
|
| 736 |
> |
It should be noted that the NPTxyz integrator is {\it not} known to |
| 737 |
> |
preserve any Hamiltonian of interest to the chemical physics |
| 738 |
> |
community. The integrator is extremely useful, however, in generating |
| 739 |
> |
initial conditions for other integration methods. It {\it is} suitable |
| 740 |
> |
for use with liquid simulations, or in cases where there is |
| 741 |
> |
orientational anisotropy in the system (i.e. in lipid bilayer |
| 742 |
> |
simulations). |
| 743 |
> |
|
| 744 |
|
\subsection{\label{Sec:zcons}Z-Constraint Method} |
| 745 |
|
|
| 746 |
< |
Based on fluctuatin-dissipation theorem,\bigskip\ force auto-correlation |
| 747 |
< |
method was developed to investigate the dynamics of ions inside the ion |
| 748 |
< |
channels.\cite{Roux91} Time-dependent friction coefficient can be calculated |
| 749 |
< |
from the deviation of the instaneous force from its mean force. |
| 746 |
> |
Based on fluctuatin-dissipation theorem,\bigskip\ force |
| 747 |
> |
auto-correlation method was developed to investigate the dynamics of |
| 748 |
> |
ions inside the ion channels.\cite{Roux91} Time-dependent friction |
| 749 |
> |
coefficient can be calculated from the deviation of the instaneous |
| 750 |
> |
force from its mean force. |
| 751 |
|
|
| 752 |
|
% |
| 753 |
|
|
