| 18 |
|
\section{\label{oopseSec:foreword}Foreword} |
| 19 |
|
|
| 20 |
|
In this chapter, I present and detail the capabilities of the open |
| 21 |
< |
source simulation package {\sc oopse}. It is important to note, that a |
| 22 |
< |
simulation package of this size and scope would not have been possible |
| 21 |
> |
source simulation program {\sc oopse}. It is important to note that a |
| 22 |
> |
simulation program of this size and scope would not have been possible |
| 23 |
|
without the collaborative efforts of my colleagues: Charles |
| 24 |
|
F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J.~Daniel |
| 25 |
|
Gezelter. Although my contributions to {\sc oopse} are major, |
| 26 |
|
consideration of my work apart from the others would not give a |
| 27 |
< |
complete description to the package's capabilities. As such, all |
| 27 |
> |
complete description to the program's capabilities. As such, all |
| 28 |
|
contributions to {\sc oopse} to date are presented in this chapter. |
| 29 |
|
|
| 30 |
|
Charles Vardeman is responsible for the parallelization of the long |
| 70 |
|
|
| 71 |
|
Despite their utility, problems with these packages arise when |
| 72 |
|
researchers try to develop techniques or energetic models that the |
| 73 |
< |
code was not originally designed to simulate. Examples of uncommonly |
| 74 |
< |
implemented techniques and energetics include; dipole-dipole |
| 75 |
< |
interactions, rigid body dynamics, and metallic embedded |
| 76 |
< |
potentials. When faced with these obstacles, a researcher must either |
| 77 |
< |
develop their own code or license and extend one of the commercial |
| 78 |
< |
packages. What we have elected to do, is develop a package of |
| 79 |
< |
simulation code capable of implementing the types of models upon which |
| 80 |
< |
our research is based. |
| 73 |
> |
code was not originally designed to simulate. Examples of techniques |
| 74 |
> |
and energetics not commonly implemented include; dipole-dipole |
| 75 |
> |
interactions, rigid body dynamics, and metallic potentials. When faced |
| 76 |
> |
with these obstacles, a researcher must either develop their own code |
| 77 |
> |
or license and extend one of the commercial packages. What we have |
| 78 |
> |
elected to do is develop a body of simulation code capable of |
| 79 |
> |
implementing the types of models upon which our research is based. |
| 80 |
|
|
| 81 |
|
In developing {\sc oopse}, we have adhered to the precepts of Open |
| 82 |
|
Source development, and are releasing our source code with a |
| 160 |
|
\begin{equation} |
| 161 |
|
\boldsymbol{\tau}_i= |
| 162 |
|
\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} |
| 163 |
< |
+ \boldsymbol{\tau}_{ia}\biggr] |
| 163 |
> |
+ \boldsymbol{\tau}_{ia}\biggr], |
| 164 |
|
\label{eq:torqueAccumulate} |
| 165 |
|
\end{equation} |
| 166 |
|
where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and |
| 172 |
|
each rigid body. In order to move between the space fixed and body |
| 173 |
|
fixed coordinate axes, parameters describing the orientation must be |
| 174 |
|
maintained for each rigid body. At a minimum, the rotation matrix |
| 175 |
< |
(\textbf{A}) can be described by the three Euler angles ($\phi, |
| 176 |
< |
\theta,$ and $\psi$), where the elements of \textbf{A} are composed of |
| 175 |
> |
($\mathsf{A}$) can be described by the three Euler angles ($\phi, |
| 176 |
> |
\theta,$ and $\psi$), where the elements of $\mathsf{A}$ are composed of |
| 177 |
|
trigonometric operations involving $\phi, \theta,$ and |
| 178 |
|
$\psi$.\cite{Goldstein01} In order to avoid numerical instabilities |
| 179 |
|
inherent in using the Euler angles, the four parameter ``quaternion'' |
| 180 |
< |
scheme is often used. The elements of \textbf{A} can be expressed as |
| 180 |
> |
scheme is often used. The elements of $\mathsf{A}$ can be expressed as |
| 181 |
|
arithmetic operations involving the four quaternions ($q_0, q_1, q_2,$ |
| 182 |
|
and $q_3$).\cite{allen87:csl} Use of quaternions also leads to |
| 183 |
|
performance enhancements, particularly for very small |
| 193 |
|
|
| 194 |
|
\begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}] |
| 195 |
|
molecule{ |
| 196 |
< |
name = "TIP3P_water"; |
| 196 |
> |
name = "TIP3P"; |
| 197 |
> |
nAtoms = 3; |
| 198 |
> |
atom[0]{ |
| 199 |
> |
type = "O_TIP3P"; |
| 200 |
> |
position( 0.0, 0.0, -0.06556 ); |
| 201 |
> |
} |
| 202 |
> |
atom[1]{ |
| 203 |
> |
type = "H_TIP3P"; |
| 204 |
> |
position( 0.0, 0.75695, 0.52032 ); |
| 205 |
> |
} |
| 206 |
> |
atom[2]{ |
| 207 |
> |
type = "H_TIP3P"; |
| 208 |
> |
position( 0.0, -0.75695, 0.52032 ); |
| 209 |
> |
} |
| 210 |
> |
|
| 211 |
|
nRigidBodies = 1; |
| 212 |
< |
rigidBody[0]{ |
| 213 |
< |
nAtoms = 3; |
| 214 |
< |
atom[0]{ |
| 202 |
< |
type = "O_TIP3P"; |
| 203 |
< |
position( 0.0, 0.0, -0.06556 ); |
| 204 |
< |
} |
| 205 |
< |
atom[1]{ |
| 206 |
< |
type = "H_TIP3P"; |
| 207 |
< |
position( 0.0, 0.75695, 0.52032 ); |
| 208 |
< |
} |
| 209 |
< |
atom[2]{ |
| 210 |
< |
type = "H_TIP3P"; |
| 211 |
< |
position( 0.0, -0.75695, 0.52032 ); |
| 212 |
< |
} |
| 213 |
< |
position( 0.0, 0.0, 0.0 ); |
| 214 |
< |
orientation( 0.0, 0.0, 1.0 ); |
| 212 |
> |
rigidBody[0]{ |
| 213 |
> |
nMembers = 3; |
| 214 |
> |
members(0, 1, 2); |
| 215 |
|
} |
| 216 |
|
} |
| 217 |
|
\end{lstlisting} |
| 227 |
|
4\epsilon_{ij} \biggl[ |
| 228 |
|
\biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{12} |
| 229 |
|
- \biggl(\frac{\sigma_{ij}}{r_{ij}}\biggr)^{6} |
| 230 |
< |
\biggr] |
| 230 |
> |
\biggr], |
| 231 |
|
\label{eq:lennardJonesPot} |
| 232 |
|
\end{equation} |
| 233 |
< |
Where $r_{ij}$ is the distance between particles $i$ and $j$, |
| 233 |
> |
where $r_{ij}$ is the distance between particles $i$ and $j$, |
| 234 |
|
$\sigma_{ij}$ scales the length of the interaction, and |
| 235 |
|
$\epsilon_{ij}$ scales the well depth of the potential. Scheme |
| 236 |
< |
\ref{sch:LJFF} gives and example \texttt{.bass} file that |
| 236 |
> |
\ref{sch:LJFF} gives an example \texttt{.bass} file that |
| 237 |
|
sets up a system of 108 Ar particles to be simulated using the |
| 238 |
|
Lennard-Jones force field. |
| 239 |
|
|
| 264 |
|
the energy value at $r_{\text{cut}}$ is subtracted from the |
| 265 |
|
potential. This causes the potential to go to zero smoothly at the |
| 266 |
|
cut-off radius, and preserves conservation of energy in integrating |
| 267 |
< |
the equations of motion. |
| 267 |
> |
the equations of motion. There still remains a discontinuity in the derivative (the forces), however, this does not significantly affect the dynamics. |
| 268 |
|
|
| 269 |
|
Interactions between dissimilar particles requires the generation of |
| 270 |
|
cross term parameters for $\sigma$ and $\epsilon$. These are |
| 271 |
|
calculated through the Lorentz-Berthelot mixing |
| 272 |
|
rules:\cite{allen87:csl} |
| 273 |
|
\begin{equation} |
| 274 |
< |
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}] |
| 274 |
> |
\sigma_{ij} = \frac{1}{2}[\sigma_{ii} + \sigma_{jj}], |
| 275 |
|
\label{eq:sigmaMix} |
| 276 |
|
\end{equation} |
| 277 |
|
and |
| 278 |
|
\begin{equation} |
| 279 |
< |
\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}} |
| 279 |
> |
\epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}. |
| 280 |
|
\label{eq:epsilonMix} |
| 281 |
|
\end{equation} |
| 282 |
|
|
| 299 |
|
include a reaction field to mimic larger range interactions. |
| 300 |
|
|
| 301 |
|
As an example, lipid head-groups in {\sc duff} are represented as |
| 302 |
< |
point dipole interaction sites. By placing a dipole at the head group |
| 303 |
< |
center of mass, our model mimics the charge separation found in common |
| 304 |
< |
phospholipids such as phosphatidylcholine.\cite{Cevc87} Additionally, |
| 305 |
< |
a large Lennard-Jones site is located at the pseudoatom's center of |
| 306 |
< |
mass. The model is illustrated by the red atom in |
| 307 |
< |
Fig.~\ref{oopseFig:lipidModel}. The water model we use to complement |
| 308 |
< |
the dipoles of the lipids is our reparameterization of the soft sticky |
| 309 |
< |
dipole (SSD) model of Ichiye |
| 302 |
> |
point dipole interaction sites. By placing a dipole at the head |
| 303 |
> |
group's center of mass, our model mimics the charge separation found |
| 304 |
> |
in common phospholipid head groups such as |
| 305 |
> |
phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones |
| 306 |
> |
site is located at the pseudoatom's center of mass. The model is |
| 307 |
> |
illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The |
| 308 |
> |
water model we use to complement the dipoles of the lipids is our |
| 309 |
> |
reparameterization of the soft sticky dipole (SSD) model of Ichiye |
| 310 |
|
\emph{et al.}\cite{liu96:new_model} |
| 311 |
|
|
| 312 |
|
\begin{figure} |
| 313 |
|
\centering |
| 314 |
< |
\includegraphics[width=\linewidth]{lipidModel.eps} |
| 314 |
> |
\includegraphics[width=\linewidth]{twoChainFig.eps} |
| 315 |
|
\caption[A representation of a lipid model in {\sc duff}]{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ % |
| 316 |
< |
is the bend angle, $\mu$ is the dipole moment of the head group, and n |
| 317 |
< |
is the chain length.} |
| 316 |
> |
is the bend angle, and $\mu$ is the dipole moment of the head group.} |
| 317 |
|
\label{oopseFig:lipidModel} |
| 318 |
|
\end{figure} |
| 319 |
|
|
| 337 |
|
integrating the equations of motion. A simulation using {\sc duff} is |
| 338 |
|
illustrated in Scheme \ref{sch:DUFF}. |
| 339 |
|
|
| 340 |
< |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}] |
| 340 |
> |
\begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion of a \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}] |
| 341 |
|
|
| 342 |
|
#include "water.mdl" |
| 343 |
|
#include "lipid.mdl" |
| 364 |
|
The total potential energy function in {\sc duff} is |
| 365 |
|
\begin{equation} |
| 366 |
|
V = \sum^{N}_{I=1} V^{I}_{\text{Internal}} |
| 367 |
< |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}} |
| 367 |
> |
+ \sum^{N-1}_{I=1} \sum_{J>I} V^{IJ}_{\text{Cross}}, |
| 368 |
|
\label{eq:totalPotential} |
| 369 |
|
\end{equation} |
| 370 |
< |
Where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
| 370 |
> |
where $V^{I}_{\text{Internal}}$ is the internal potential of molecule $I$: |
| 371 |
|
\begin{equation} |
| 372 |
|
V^{I}_{\text{Internal}} = |
| 373 |
|
\sum_{\theta_{ijk} \in I} V_{\text{bend}}(\theta_{ijk}) |
| 375 |
|
+ \sum_{i \in I} \sum_{(j>i+4) \in I} |
| 376 |
|
\biggl[ V_{\text{LJ}}(r_{ij}) + V_{\text{dipole}} |
| 377 |
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 378 |
< |
\biggr] |
| 378 |
> |
\biggr]. |
| 379 |
|
\label{eq:internalPotential} |
| 380 |
|
\end{equation} |
| 381 |
|
Here $V_{\text{bend}}$ is the bend potential for all 1, 3 bonded pairs |
| 382 |
|
within the molecule $I$, and $V_{\text{torsion}}$ is the torsion potential |
| 383 |
|
for all 1, 4 bonded pairs. The pairwise portions of the internal |
| 384 |
< |
potential are excluded for pairs that are closer than three bonds, |
| 386 |
< |
i.e.~atom pairs farther away than a torsion are included in the |
| 387 |
< |
pair-wise loop. |
| 384 |
> |
potential are excluded for atom pairs that are involved in the same bond, bend, or torsion. All other atom pairs within the molecule are subject to the LJ pair potential. |
| 385 |
|
|
| 386 |
|
|
| 387 |
|
The bend potential of a molecule is represented by the following function: |
| 388 |
|
\begin{equation} |
| 389 |
< |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2 \label{eq:bendPot} |
| 389 |
> |
V_{\text{bend}}(\theta_{ijk}) = k_{\theta}( \theta_{ijk} - \theta_0 )^2, \label{eq:bendPot} |
| 390 |
|
\end{equation} |
| 391 |
< |
Where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
| 391 |
> |
where $\theta_{ijk}$ is the angle defined by atoms $i$, $j$, and $k$ |
| 392 |
|
(see Fig.~\ref{oopseFig:lipidModel}), $\theta_0$ is the equilibrium |
| 393 |
|
bond angle, and $k_{\theta}$ is the force constant which determines the |
| 394 |
|
strength of the harmonic bend. The parameters for $k_{\theta}$ and |
| 399 |
|
\begin{equation} |
| 400 |
|
V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi] |
| 401 |
|
+ c_2[1 + \cos(2\phi)] |
| 402 |
< |
+ c_3[1 + \cos(3\phi)] |
| 402 |
> |
+ c_3[1 + \cos(3\phi)], |
| 403 |
|
\label{eq:origTorsionPot} |
| 404 |
|
\end{equation} |
| 405 |
< |
Where: |
| 405 |
> |
where: |
| 406 |
|
\begin{equation} |
| 407 |
|
\cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot |
| 408 |
< |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}) |
| 408 |
> |
(\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl}). |
| 409 |
|
\label{eq:torsPhi} |
| 410 |
|
\end{equation} |
| 411 |
|
Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond |
| 415 |
|
a power series of the form: |
| 416 |
|
\begin{equation} |
| 417 |
|
V_{\text{torsion}}(\phi) = |
| 418 |
< |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0 |
| 418 |
> |
k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0, |
| 419 |
|
\label{eq:torsionPot} |
| 420 |
|
\end{equation} |
| 421 |
< |
Where: |
| 421 |
> |
where: |
| 422 |
|
\begin{align*} |
| 423 |
< |
k_0 &= c_1 + c_3 \\ |
| 424 |
< |
k_1 &= c_1 - 3c_3 \\ |
| 425 |
< |
k_2 &= 2 c_2 \\ |
| 426 |
< |
k_3 &= 4c_3 |
| 423 |
> |
k_0 &= c_1 + c_3, \\ |
| 424 |
> |
k_1 &= c_1 - 3c_3, \\ |
| 425 |
> |
k_2 &= 2 c_2, \\ |
| 426 |
> |
k_3 &= 4c_3. |
| 427 |
|
\end{align*} |
| 428 |
|
By recasting the potential as a power series, repeated trigonometric |
| 429 |
|
evaluations are avoided during the calculation of the potential energy. |
| 438 |
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 439 |
|
+ V_{\text{sticky}} |
| 440 |
|
(\mathbf{r}_{ij},\boldsymbol{\Omega}_{i},\boldsymbol{\Omega}_{j}) |
| 441 |
< |
\biggr] |
| 441 |
> |
\biggr], |
| 442 |
|
\label{eq:crossPotentail} |
| 443 |
|
\end{equation} |
| 444 |
< |
Where $V_{\text{LJ}}$ is the Lennard Jones potential, |
| 444 |
> |
where $V_{\text{LJ}}$ is the Lennard Jones potential, |
| 445 |
|
$V_{\text{dipole}}$ is the dipole dipole potential, and |
| 446 |
|
$V_{\text{sticky}}$ is the sticky potential defined by the SSD model |
| 447 |
|
(Sec.~\ref{oopseSec:SSD}). Note that not all atom types include all |
| 454 |
|
\boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j} |
| 455 |
|
- |
| 456 |
|
3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) % |
| 457 |
< |
(\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr] |
| 457 |
> |
(\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr]. |
| 458 |
|
\label{eq:dipolePot} |
| 459 |
|
\end{equation} |
| 460 |
|
Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing |
| 492 |
|
|
| 493 |
|
In the interest of computational efficiency, the default solvent used |
| 494 |
|
by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water |
| 495 |
< |
model.\cite{Gezelter04} The original SSD was developed by Ichiye |
| 495 |
> |
model.\cite{fennell04} The original SSD was developed by Ichiye |
| 496 |
|
\emph{et al.}\cite{liu96:new_model} as a modified form of the hard-sphere |
| 497 |
|
water model proposed by Bratko, Blum, and |
| 498 |
|
Luzar.\cite{Bratko85,Bratko95} It consists of a single point dipole |
| 566 |
|
|
| 567 |
|
Recent constant pressure simulations revealed issues in the original |
| 568 |
|
SSD model that led to lower than expected densities at all target |
| 569 |
< |
pressures.\cite{Ichiye03,Gezelter04} The default model in {\sc oopse} |
| 569 |
> |
pressures.\cite{Ichiye03,fennell04} The default model in {\sc oopse} |
| 570 |
|
is therefore SSD/E, a density corrected derivative of SSD that |
| 571 |
|
exhibits improved liquid structure and transport behavior. If the use |
| 572 |
|
of a reaction field long-range interaction correction is desired, it |
| 573 |
|
is recommended that the parameters be modified to those of the SSD/RF |
| 574 |
< |
model. Solvent parameters can be easily modified in an accompanying |
| 574 |
> |
model (an SSD variant parameterized for reaction field). Solvent parameters can be easily modified in an accompanying |
| 575 |
|
\texttt{.bass} file as illustrated in the scheme below. A table of the |
| 576 |
|
parameter values and the drawbacks and benefits of the different |
| 577 |
|
density corrected SSD models can be found in |
| 578 |
< |
reference~\cite{Gezelter04}. |
| 578 |
> |
reference~\cite{fennell04}. |
| 579 |
|
|
| 580 |
< |
\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]An example file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
| 580 |
> |
\begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]A portion of a \texttt{.bass} file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
| 581 |
|
|
| 582 |
|
#include "water.mdl" |
| 583 |
|
|
| 625 |
|
The {\sc eam} potential has the form: |
| 626 |
|
\begin{eqnarray} |
| 627 |
|
V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i} |
| 628 |
< |
\phi_{ij}({\bf r}_{ij}) \\ |
| 629 |
< |
\rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}) |
| 628 |
> |
\phi_{ij}({\bf r}_{ij}), \\ |
| 629 |
> |
\rho_{i} & = & \sum_{j \neq i} f_{j}({\bf r}_{ij}), |
| 630 |
|
\end{eqnarray} |
| 631 |
|
where $F_{i} $ is the embedding function that equates the energy |
| 632 |
|
required to embed a positively-charged core ion $i$ into a linear |
| 641 |
|
surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$ |
| 642 |
|
interactions. Foiles \emph{et al}.~fit {\sc eam} potentials for the fcc |
| 643 |
|
metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86} |
| 644 |
< |
These fits, are included in {\sc oopse}. |
| 644 |
> |
These fits are included in {\sc oopse}. |
| 645 |
|
|
| 646 |
|
\subsection{\label{oopseSec:pbc}Periodic Boundary Conditions} |
| 647 |
|
|
| 659 |
|
use a $3 \times 3$ matrix, $\mathsf{H}$, to describe the shape and |
| 660 |
|
size of the simulation box. $\mathsf{H}$ is defined: |
| 661 |
|
\begin{equation} |
| 662 |
< |
\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ) |
| 662 |
> |
\mathsf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z ), |
| 663 |
|
\end{equation} |
| 664 |
< |
Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the |
| 664 |
> |
where $\mathbf{h}_{\alpha}$ is the column vector of the $\alpha$ axis of the |
| 665 |
|
box. During the course of the simulation both the size and shape of |
| 666 |
|
the box can be changed to allow volume fluctuations when constraining |
| 667 |
|
the pressure. |
| 669 |
|
A real space vector, $\mathbf{r}$ can be transformed in to a box space |
| 670 |
|
vector, $\mathbf{s}$, and back through the following transformations: |
| 671 |
|
\begin{align} |
| 672 |
< |
\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r} \\ |
| 673 |
< |
\mathbf{r} &= \mathsf{H} \mathbf{s} |
| 672 |
> |
\mathbf{s} &= \mathsf{H}^{-1} \mathbf{r}, \\ |
| 673 |
> |
\mathbf{r} &= \mathsf{H} \mathbf{s}. |
| 674 |
|
\end{align} |
| 675 |
|
The vector $\mathbf{s}$ is now a vector expressed as the number of box |
| 676 |
|
lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$ |
| 677 |
|
directions. To find the minimum image of a vector $\mathbf{r}$, we |
| 678 |
|
first convert it to its corresponding vector in box space, and then, |
| 679 |
< |
cast each element to lie on the in the range $[-0.5,0.5]$: |
| 679 |
> |
cast each element to lie in the range $[-0.5,0.5]$: |
| 680 |
|
\begin{equation} |
| 681 |
< |
s_{i}^{\prime}=s_{i}-\roundme(s_{i}) |
| 681 |
> |
s_{i}^{\prime}=s_{i}-\roundme(s_{i}), |
| 682 |
|
\end{equation} |
| 683 |
< |
Where $s_i$ is the $i$th element of $\mathbf{s}$, and |
| 684 |
< |
$\roundme(s_i)$is given by |
| 683 |
> |
where $s_i$ is the $i$th element of $\mathbf{s}$, and |
| 684 |
> |
$\roundme(s_i)$ is given by |
| 685 |
|
\begin{equation} |
| 686 |
|
\roundme(x) = |
| 687 |
|
\begin{cases} |
| 688 |
< |
\lfloor x+0.5 \rfloor & \text{if $x \ge 0$} \\ |
| 689 |
< |
\lceil x-0.5 \rceil & \text{if $x < 0$ } |
| 688 |
> |
\lfloor x+0.5 \rfloor & \text{if $x \ge 0$,} \\ |
| 689 |
> |
\lceil x-0.5 \rceil & \text{if $x < 0$.} |
| 690 |
|
\end{cases} |
| 691 |
|
\end{equation} |
| 692 |
|
Here $\lfloor x \rfloor$ is the floor operator, and gives the largest |
| 698 |
|
Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by |
| 699 |
|
transforming back to real space, |
| 700 |
|
\begin{equation} |
| 701 |
< |
\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}% |
| 701 |
> |
\mathbf{r}^{\prime}=\mathsf{H}^{-1}\mathbf{s}^{\prime}.% |
| 702 |
|
\end{equation} |
| 703 |
|
In this way, particles are allowed to diffuse freely in $\mathbf{r}$, |
| 704 |
|
but their minimum images, $\mathbf{r}^{\prime}$ are used to compute |
| 808 |
|
entities are written out using quanternions, to save space in the |
| 809 |
|
output files. |
| 810 |
|
|
| 811 |
< |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well.. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular velocities.},label=sch:dumpFormat] |
| 811 |
> |
\begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular velocities.},label=sch:dumpFormat] |
| 812 |
|
|
| 813 |
|
nAtoms |
| 814 |
|
time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz; |
| 864 |
|
|
| 865 |
|
\section{\label{oopseSec:mechanics}Mechanics} |
| 866 |
|
|
| 870 |
– |
|
| 871 |
– |
\section{\label{sec:mechanics}Mechanics} |
| 872 |
– |
|
| 867 |
|
\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the |
| 868 |
|
DLM method} |
| 869 |
|
|
| 873 |
|
(DLM).\cite{Dullweber1997} When there are no directional atoms or |
| 874 |
|
rigid bodies present in the simulation, this integrator becomes the |
| 875 |
|
standard velocity-Verlet integrator which is known to sample the |
| 876 |
< |
microcanonical (NVE) ensemble.\cite{} |
| 876 |
> |
microcanonical (NVE) ensemble.\cite{Frenkel1996} |
| 877 |
|
|
| 878 |
|
Previous integration methods for orientational motion have problems |
| 879 |
|
that are avoided in the DLM method. Direct propagation of the Euler |
| 906 |
|
H = \sum_{i} \left( \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 907 |
|
\frac{1}{2} {\bf j}_i^T \cdot \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot |
| 908 |
|
{\bf j}_i \right) + |
| 909 |
< |
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right) |
| 909 |
> |
V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right), |
| 910 |
|
\end{equation} |
| 911 |
|
where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector |
| 912 |
< |
and velocity of the center of mass of particle $i$, and ${\bf j}_i$ |
| 913 |
< |
and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
| 914 |
< |
momentum and moment of inertia tensor, respectively. $\mathsf{A}_i$ |
| 912 |
> |
and velocity of the center of mass of particle $i$, and ${\bf j}_i$, |
| 913 |
> |
$\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
| 914 |
> |
momentum and moment of inertia tensor respectively, and the |
| 915 |
> |
superscript $T$ denotes the transpose of the vector. $\mathsf{A}_i$ |
| 916 |
|
is the $3 \times 3$ rotation matrix describing the instantaneous |
| 917 |
|
orientation of the particle. $V$ is the potential energy function |
| 918 |
|
which may depend on both the positions $\left\{{\bf r}\right\}$ and |
| 920 |
|
equations of motion for the particle centers of mass are derived from |
| 921 |
|
Hamilton's equations and are quite simple, |
| 922 |
|
\begin{eqnarray} |
| 923 |
< |
\dot{{\bf r}} & = & {\bf v} \\ |
| 924 |
< |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} |
| 923 |
> |
\dot{{\bf r}} & = & {\bf v}, \\ |
| 924 |
> |
\dot{{\bf v}} & = & \frac{{\bf f}}{m}, |
| 925 |
|
\end{eqnarray} |
| 926 |
|
where ${\bf f}$ is the instantaneous force on the center of mass |
| 927 |
|
of the particle, |
| 933 |
|
The equations of motion for the orientational degrees of freedom are |
| 934 |
|
\begin{eqnarray} |
| 935 |
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 936 |
< |
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\ |
| 936 |
> |
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right),\\ |
| 937 |
|
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} |
| 938 |
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 939 |
< |
V}{\partial \mathsf{A}} \right) |
| 939 |
> |
V}{\partial \mathsf{A}} \right). |
| 940 |
|
\end{eqnarray} |
| 941 |
|
In these equations of motion, the $\mbox{skew}$ matrix of a vector |
| 942 |
|
${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined: |
| 947 |
|
-v_3 & 0 & v_1 \\ |
| 948 |
|
v_2 & -v_1 & 0 |
| 949 |
|
\end{array} |
| 950 |
< |
\right) |
| 950 |
> |
\right). |
| 951 |
|
\end{equation} |
| 952 |
|
The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$ |
| 953 |
|
rotation matrix to a vector of orientations by first computing the |
| 956 |
|
$\mbox{skew}$ function above: |
| 957 |
|
\begin{equation} |
| 958 |
|
\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A} |
| 959 |
< |
- \mathsf{A}^{T} \right) |
| 959 |
> |
- \mathsf{A}^{T} \right). |
| 960 |
|
\end{equation} |
| 961 |
|
Written this way, the $\mbox{rot}$ operation creates a set of |
| 962 |
|
conjugate angle coordinates to the body-fixed angular momenta |
| 964 |
|
is equivalent to the more familiar body-fixed forms, |
| 965 |
|
\begin{eqnarray} |
| 966 |
|
\dot{j_{x}} & = & \tau^b_x(t) + |
| 967 |
< |
\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z \\ |
| 967 |
> |
\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z, \\ |
| 968 |
|
\dot{j_{y}} & = & \tau^b_y(t) + |
| 969 |
< |
\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x \\ |
| 969 |
> |
\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x,\\ |
| 970 |
|
\dot{j_{z}} & = & \tau^b_z(t) + |
| 971 |
< |
\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y |
| 971 |
> |
\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y, |
| 972 |
|
\end{eqnarray} |
| 973 |
|
which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are |
| 974 |
|
most easily derived in the space-fixed frame, |
| 975 |
|
\begin{equation} |
| 976 |
< |
{\bf \tau}^b(t) = \mathsf{A}(t) \cdot {\bf \tau}^s(t) |
| 976 |
> |
{\bf \tau}^b(t) = \mathsf{A}(t) \cdot {\bf \tau}^s(t), |
| 977 |
|
\end{equation} |
| 978 |
|
where the torques are either derived from the forces on the |
| 979 |
|
constituent atoms of the rigid body, or for directional atoms, |
| 993 |
|
{\it symplectic}), |
| 994 |
|
\item the integrator is time-{\it reversible}, making it suitable for Hybrid |
| 995 |
|
Monte Carlo applications, and |
| 996 |
< |
\item the error for a single time step is of order $O\left(h^3\right)$ |
| 996 |
> |
\item the error for a single time step is of order $\mathcal{O}\left(h^4\right)$ |
| 997 |
|
for timesteps of length $h$. |
| 998 |
|
\end{enumerate} |
| 999 |
|
|
| 1000 |
|
The integration of the equations of motion is carried out in a |
| 1001 |
< |
velocity-Verlet style 2-part algorithm: |
| 1001 |
> |
velocity-Verlet style 2-part algorithm, where $h= \delta t$: |
| 1002 |
|
|
| 1003 |
|
{\tt moveA:} |
| 1004 |
< |
\begin{eqnarray} |
| 1005 |
< |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 1006 |
< |
v}(t) + \frac{\delta t}{2} \left( {\bf f}(t) / m \right) \\ |
| 1007 |
< |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf |
| 1008 |
< |
v}\left(t + \delta t / 2 \right) \\ |
| 1009 |
< |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 1010 |
< |
j}(t) + \frac{\delta t}{2} {\bf \tau}^b(t) \\ |
| 1011 |
< |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left( \delta t |
| 1012 |
< |
{\bf j}(t + \delta t / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} |
| 1013 |
< |
\right) |
| 1014 |
< |
\end{eqnarray} |
| 1004 |
> |
\begin{align*} |
| 1005 |
> |
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1006 |
> |
+ \frac{h}{2} \left( {\bf f}(t) / m \right), \\ |
| 1007 |
> |
% |
| 1008 |
> |
{\bf r}(t + h) &\leftarrow {\bf r}(t) |
| 1009 |
> |
+ h {\bf v}\left(t + h / 2 \right), \\ |
| 1010 |
> |
% |
| 1011 |
> |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1012 |
> |
+ \frac{h}{2} {\bf \tau}^b(t), \\ |
| 1013 |
> |
% |
| 1014 |
> |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} |
| 1015 |
> |
(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right). |
| 1016 |
> |
\end{align*} |
| 1017 |
|
|
| 1018 |
< |
In this context, the $\mathrm{rot}$ function is the reversible product |
| 1018 |
> |
In this context, the $\mathrm{rotate}$ function is the reversible product |
| 1019 |
|
of the three body-fixed rotations, |
| 1020 |
|
\begin{equation} |
| 1021 |
< |
\mathrm{rot}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
| 1021 |
> |
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot |
| 1022 |
|
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y / |
| 1023 |
< |
2) \cdot \mathsf{G}_x(a_x /2) |
| 1023 |
> |
2) \cdot \mathsf{G}_x(a_x /2), |
| 1024 |
|
\end{equation} |
| 1025 |
|
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates |
| 1026 |
|
both the rotation matrix ($\mathsf{A}$) and the body-fixed angular |
| 1029 |
|
\begin{equation} |
| 1030 |
|
\mathsf{G}_\alpha( \theta ) = \left\{ |
| 1031 |
|
\begin{array}{lcl} |
| 1032 |
< |
\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T \\ |
| 1033 |
< |
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0) |
| 1032 |
> |
\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ |
| 1033 |
> |
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0). |
| 1034 |
|
\end{array} |
| 1035 |
|
\right. |
| 1036 |
|
\end{equation} |
| 1055 |
|
torques are calculated at the new positions and orientations |
| 1056 |
|
|
| 1057 |
|
{\tt doForces:} |
| 1058 |
< |
\begin{eqnarray} |
| 1059 |
< |
{\bf f}(t + \delta t) & \leftarrow & - \left(\frac{\partial V}{\partial {\bf |
| 1060 |
< |
r}}\right)_{{\bf r}(t + \delta t)} \\ |
| 1061 |
< |
{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t) |
| 1062 |
< |
\times \frac{\partial V}{\partial {\bf u}} \\ |
| 1063 |
< |
{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t) |
| 1064 |
< |
\cdot {\bf \tau}^s(t + \delta t) |
| 1065 |
< |
\end{eqnarray} |
| 1058 |
> |
\begin{align*} |
| 1059 |
> |
{\bf f}(t + h) &\leftarrow |
| 1060 |
> |
- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\ |
| 1061 |
> |
% |
| 1062 |
> |
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) |
| 1063 |
> |
\times \frac{\partial V}{\partial {\bf u}}, \\ |
| 1064 |
> |
% |
| 1065 |
> |
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h) |
| 1066 |
> |
\cdot {\bf \tau}^s(t + h). |
| 1067 |
> |
\end{align*} |
| 1068 |
|
|
| 1069 |
|
{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix |
| 1070 |
|
$\mathsf{A}$ is calculated in {\tt moveA}. Once the forces and |
| 1072 |
|
advanced to the same time value. |
| 1073 |
|
|
| 1074 |
|
{\tt moveB:} |
| 1075 |
< |
\begin{eqnarray} |
| 1076 |
< |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1077 |
< |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 1078 |
< |
{\bf f}(t + \delta t) / m \right) \\ |
| 1079 |
< |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1080 |
< |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} {\bf |
| 1081 |
< |
\tau}^b(t + \delta t) |
| 1083 |
< |
\end{eqnarray} |
| 1075 |
> |
\begin{align*} |
| 1076 |
> |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t + h / 2 \right) |
| 1077 |
> |
+ \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\ |
| 1078 |
> |
% |
| 1079 |
> |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t + h / 2 \right) |
| 1080 |
> |
+ \frac{h}{2} {\bf \tau}^b(t + h) . |
| 1081 |
> |
\end{align*} |
| 1082 |
|
|
| 1083 |
|
The matrix rotations used in the DLM method end up being more costly |
| 1084 |
|
computationally than the simpler arithmetic quaternion |
| 1086 |
|
shows an average 7\% increase in computation time using the DLM method |
| 1087 |
|
in place of quaternions. This cost is more than justified when |
| 1088 |
|
comparing the energy conservation of the two methods as illustrated in |
| 1089 |
< |
figure \ref{timestep}. |
| 1089 |
> |
Fig.~\ref{timestep}. |
| 1090 |
|
|
| 1091 |
|
\begin{figure} |
| 1092 |
|
\centering |
| 1100 |
|
\label{timestep} |
| 1101 |
|
\end{figure} |
| 1102 |
|
|
| 1103 |
< |
In figure \ref{timestep}, the resulting energy drift at various time |
| 1103 |
> |
In Fig.~\ref{timestep}, the resulting energy drift at various time |
| 1104 |
|
steps for both the DLM and quaternion integration schemes is |
| 1105 |
|
compared. All of the 1000 molecule water simulations started with the |
| 1106 |
|
same configuration, and the only difference was the method for |
| 1120 |
|
\begin{tabular}{llll} |
| 1121 |
|
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 1122 |
|
default value} \\ |
| 1123 |
< |
$\delta t$ & {\tt dt = 2.0;} & fs & none |
| 1123 |
> |
$h$ & {\tt dt = 2.0;} & fs & none |
| 1124 |
|
\end{tabular} |
| 1125 |
|
\end{center} |
| 1126 |
|
|
| 1134 |
|
\begin{center} |
| 1135 |
|
\begin{tabular}{lll} |
| 1136 |
|
{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\ |
| 1137 |
< |
NVE & microcanonical & {\tt ensemble = ``NVE''; } \\ |
| 1138 |
< |
NVT & canonical & {\tt ensemble = ``NVT''; } \\ |
| 1139 |
< |
NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt |
| 1140 |
< |
ensemble = ``NPTi'';} \\ |
| 1141 |
< |
NPTf & isobaric-isothermal (with changes to box shape) & {\tt |
| 1142 |
< |
ensemble = ``NPTf'';} \\ |
| 1143 |
< |
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = |
| 1144 |
< |
``NPTxyz'';} \\ |
| 1147 |
< |
& (with separate barostats on each box dimension) & |
| 1137 |
> |
NVE & microcanonical & {\tt ensemble = NVE; } \\ |
| 1138 |
> |
NVT & canonical & {\tt ensemble = NVT; } \\ |
| 1139 |
> |
NPTi & isobaric-isothermal & {\tt ensemble = NPTi;} \\ |
| 1140 |
> |
& (with isotropic volume changes) & \\ |
| 1141 |
> |
NPTf & isobaric-isothermal & {\tt ensemble = NPTf;} \\ |
| 1142 |
> |
& (with changes to box shape) & \\ |
| 1143 |
> |
NPTxyz & approximate isobaric-isothermal & {\tt ensemble = NPTxyz;} \\ |
| 1144 |
> |
& (with separate barostats on each box dimension) & \\ |
| 1145 |
|
\end{tabular} |
| 1146 |
|
\end{center} |
| 1147 |
|
|
| 1148 |
< |
The relatively well-known Nos\'e-Hoover thermostat is implemented in |
| 1149 |
< |
{\sc oopse}'s NVT integrator. This method couples an extra degree of |
| 1150 |
< |
freedom (the thermostat) to the kinetic energy of the system, and has |
| 1151 |
< |
been shown to sample the canonical distribution in the system degrees |
| 1152 |
< |
of freedom while conserving a quantity that is, to within a constant, |
| 1153 |
< |
the Helmholtz free energy. |
| 1148 |
> |
The relatively well-known Nos\'e-Hoover thermostat\cite{Hoover85} is |
| 1149 |
> |
implemented in {\sc oopse}'s NVT integrator. This method couples an |
| 1150 |
> |
extra degree of freedom (the thermostat) to the kinetic energy of the |
| 1151 |
> |
system, and has been shown to sample the canonical distribution in the |
| 1152 |
> |
system degrees of freedom while conserving a quantity that is, to |
| 1153 |
> |
within a constant, the Helmholtz free energy.\cite{melchionna93} |
| 1154 |
|
|
| 1155 |
|
NPT algorithms attempt to maintain constant pressure in the system by |
| 1156 |
|
coupling the volume of the system to a barostat. {\sc oopse} contains |
| 1157 |
|
three different constant pressure algorithms. The first two, NPTi and |
| 1158 |
|
NPTf have been shown to conserve a quantity that is, to within a |
| 1159 |
< |
constant, the Gibbs free energy. The Melchionna modification to the |
| 1160 |
< |
Hoover barostat is implemented in both NPTi and NPTf. NPTi allows |
| 1161 |
< |
only isotropic changes in the simulation box, while box {\it shape} |
| 1162 |
< |
variations are allowed in NPTf. The NPTxyz integrator has {\it not} |
| 1163 |
< |
been shown to sample from the isobaric-isothermal ensemble. It is |
| 1164 |
< |
useful, however, in that it maintains orthogonality for the axes of |
| 1165 |
< |
the simulation box while attempting to equalize pressure along the |
| 1166 |
< |
three perpendicular directions in the box. |
| 1159 |
> |
constant, the Gibbs free energy.\cite{melchionna93} The Melchionna |
| 1160 |
> |
modification to the Hoover barostat is implemented in both NPTi and |
| 1161 |
> |
NPTf. NPTi allows only isotropic changes in the simulation box, while |
| 1162 |
> |
box {\it shape} variations are allowed in NPTf. The NPTxyz integrator |
| 1163 |
> |
has {\it not} been shown to sample from the isobaric-isothermal |
| 1164 |
> |
ensemble. It is useful, however, in that it maintains orthogonality |
| 1165 |
> |
for the axes of the simulation box while attempting to equalize |
| 1166 |
> |
pressure along the three perpendicular directions in the box. |
| 1167 |
|
|
| 1168 |
|
Each of the extended system integrators requires additional keywords |
| 1169 |
|
to set target values for the thermodynamic state variables that are |
| 1171 |
|
characteristic decay times for the dynamics of the extended |
| 1172 |
|
variables. |
| 1173 |
|
|
| 1174 |
+ |
\begin{center} |
| 1175 |
|
\begin{tabular}{llll} |
| 1176 |
|
{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf |
| 1177 |
|
default value} \\ |
| 1180 |
|
$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\ |
| 1181 |
|
$\tau_B$ & {\tt tauBarostat = 5e3;} & fs & none \\ |
| 1182 |
|
& {\tt resetTime = 200;} & fs & none \\ |
| 1183 |
< |
& {\tt useInitialExtendedSystemState = ``true'';} & logical & |
| 1184 |
< |
false |
| 1183 |
> |
& {\tt useInitialExtendedSystemState = true;} & logical & |
| 1184 |
> |
true |
| 1185 |
|
\end{tabular} |
| 1186 |
+ |
\end{center} |
| 1187 |
|
|
| 1188 |
|
Two additional keywords can be used to either clear the extended |
| 1189 |
|
system variables periodically ({\tt resetTime}), or to maintain the |
| 1191 |
|
useInitialExtendedSystemState}). More details on these variables |
| 1192 |
|
and their use in the integrators follows below. |
| 1193 |
|
|
| 1194 |
< |
\subsubsection{\label{oopseSec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting} |
| 1194 |
> |
\subsection{\label{oopseSec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting} |
| 1195 |
|
|
| 1196 |
|
The Nos\'e-Hoover equations of motion are given by\cite{Hoover85} |
| 1197 |
|
\begin{eqnarray} |
| 1198 |
< |
\dot{{\bf r}} & = & {\bf v} \\ |
| 1199 |
< |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\ |
| 1198 |
> |
\dot{{\bf r}} & = & {\bf v}, \\ |
| 1199 |
> |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} ,\\ |
| 1200 |
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1201 |
< |
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\ |
| 1201 |
> |
\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right), \\ |
| 1202 |
|
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1} |
| 1203 |
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1204 |
< |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} |
| 1204 |
> |
V}{\partial \mathsf{A}} \right) - \chi {\bf j}. |
| 1205 |
|
\label{eq:nosehoovereom} |
| 1206 |
|
\end{eqnarray} |
| 1207 |
|
|
| 1219 |
|
\end{equation} |
| 1220 |
|
Here, $f$ is the total number of degrees of freedom in the system, |
| 1221 |
|
\begin{equation} |
| 1222 |
< |
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}} |
| 1222 |
> |
f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}}, |
| 1223 |
|
\end{equation} |
| 1224 |
|
and $K$ is the total kinetic energy, |
| 1225 |
|
\begin{equation} |
| 1226 |
|
K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i + |
| 1227 |
|
\sum_{i=1}^{N_{\mathrm{orient}}} \frac{1}{2} {\bf j}_i^T \cdot |
| 1228 |
< |
\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i |
| 1228 |
> |
\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i. |
| 1229 |
|
\end{equation} |
| 1230 |
|
|
| 1231 |
|
In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for |
| 1238 |
|
part algorithm: |
| 1239 |
|
|
| 1240 |
|
{\tt moveA:} |
| 1241 |
< |
\begin{eqnarray} |
| 1242 |
< |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 1243 |
< |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 1244 |
< |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1245 |
< |
\chi(t)\right) \\ |
| 1246 |
< |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf |
| 1247 |
< |
v}\left(t + \delta t / 2 \right) \\ |
| 1248 |
< |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 1249 |
< |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1250 |
< |
\chi(t) \right) \\ |
| 1251 |
< |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 1252 |
< |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\ |
| 1253 |
< |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 1254 |
< |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1255 |
< |
\right) |
| 1256 |
< |
\end{eqnarray} |
| 1241 |
> |
\begin{align*} |
| 1242 |
> |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1243 |
> |
% |
| 1244 |
> |
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1245 |
> |
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1246 |
> |
\chi(t)\right), \\ |
| 1247 |
> |
% |
| 1248 |
> |
{\bf r}(t + h) &\leftarrow {\bf r}(t) |
| 1249 |
> |
+ h {\bf v}\left(t + h / 2 \right) ,\\ |
| 1250 |
> |
% |
| 1251 |
> |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1252 |
> |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1253 |
> |
\chi(t) \right) ,\\ |
| 1254 |
> |
% |
| 1255 |
> |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate} |
| 1256 |
> |
\left(h * {\bf j}(t + h / 2) |
| 1257 |
> |
\overleftrightarrow{\mathsf{I}}^{-1} \right) ,\\ |
| 1258 |
> |
% |
| 1259 |
> |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) |
| 1260 |
> |
+ \frac{h}{2 \tau_T^2} \left( \frac{T(t)} |
| 1261 |
> |
{T_{\mathrm{target}}} - 1 \right) . |
| 1262 |
> |
\end{align*} |
| 1263 |
|
|
| 1264 |
< |
Here $\mathrm{rot}(\delta t * {\bf j} |
| 1264 |
> |
Here $\mathrm{rotate}(h * {\bf j} |
| 1265 |
|
\overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Trotter |
| 1266 |
|
factorization of the three rotation operations that was discussed in |
| 1267 |
|
the section on the DLM integrator. Note that this operation modifies |
| 1277 |
|
advanced to the same time value. |
| 1278 |
|
|
| 1279 |
|
{\tt moveB:} |
| 1280 |
< |
\begin{eqnarray} |
| 1281 |
< |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 1282 |
< |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 1283 |
< |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 1284 |
< |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 1285 |
< |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 1286 |
< |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1287 |
< |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 1288 |
< |
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t) |
| 1289 |
< |
\chi(t \delta t)\right) \\ |
| 1290 |
< |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1291 |
< |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 1292 |
< |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 1293 |
< |
\chi(t + \delta t) \right) |
| 1294 |
< |
\end{eqnarray} |
| 1280 |
> |
\begin{align*} |
| 1281 |
> |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1282 |
> |
\left\{{\bf j}(t + h)\right\}, \\ |
| 1283 |
> |
% |
| 1284 |
> |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1285 |
> |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} |
| 1286 |
> |
{T_{\mathrm{target}}} - 1 \right), \\ |
| 1287 |
> |
% |
| 1288 |
> |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t |
| 1289 |
> |
+ h / 2 \right) + \frac{h}{2} \left( |
| 1290 |
> |
\frac{{\bf f}(t + h)}{m} - {\bf v}(t + h) |
| 1291 |
> |
\chi(t h)\right) ,\\ |
| 1292 |
> |
% |
| 1293 |
> |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
| 1294 |
> |
+ h / 2 \right) + \frac{h}{2} |
| 1295 |
> |
\left( {\bf \tau}^b(t + h) - {\bf j}(t + h) |
| 1296 |
> |
\chi(t + h) \right) . |
| 1297 |
> |
\end{align*} |
| 1298 |
|
|
| 1299 |
< |
Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required |
| 1300 |
< |
to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, they |
| 1301 |
< |
indirectly depend on their own values at time $t + \delta t$. {\tt |
| 1302 |
< |
moveB} is therefore done in an iterative fashion until $\chi(t + |
| 1303 |
< |
\delta t)$ becomes self-consistent. The relative tolerance for the |
| 1304 |
< |
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but |
| 1305 |
< |
{\sc oopse} will terminate the iteration after 4 loops even if the |
| 1298 |
< |
consistency check has not been satisfied. |
| 1299 |
> |
Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to caclculate |
| 1300 |
> |
$T(t + h)$ as well as $\chi(t + h)$, they indirectly depend on their |
| 1301 |
> |
own values at time $t + h$. {\tt moveB} is therefore done in an |
| 1302 |
> |
iterative fashion until $\chi(t + h)$ becomes self-consistent. The |
| 1303 |
> |
relative tolerance for the self-consistency check defaults to a value |
| 1304 |
> |
of $\mbox{10}^{-6}$, but {\sc oopse} will terminate the iteration |
| 1305 |
> |
after 4 loops even if the consistency check has not been satisfied. |
| 1306 |
|
|
| 1307 |
|
The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the |
| 1308 |
|
extended system that is, to within a constant, identical to the |
| 1309 |
< |
Helmholtz free energy, |
| 1309 |
> |
Helmholtz free energy,\cite{melchionna93} |
| 1310 |
|
\begin{equation} |
| 1311 |
|
H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left( |
| 1312 |
|
\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime |
| 1313 |
< |
\right) |
| 1313 |
> |
\right). |
| 1314 |
|
\end{equation} |
| 1315 |
< |
Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation |
| 1315 |
> |
Poor choices of $h$ or $\tau_T$ can result in non-conservation |
| 1316 |
|
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the |
| 1317 |
|
last column of the {\tt .stat} file to allow checks on the quality of |
| 1318 |
|
the integration. |
| 1321 |
|
{\tt moveB} portions of the algorithm. Details on the constraint |
| 1322 |
|
algorithms are given in section \ref{oopseSec:rattle}. |
| 1323 |
|
|
| 1324 |
< |
\subsubsection{\label{sec:NPTi}Constant-pressure integration with |
| 1324 |
> |
\subsection{\label{sec:NPTi}Constant-pressure integration with |
| 1325 |
|
isotropic box deformations (NPTi)} |
| 1326 |
|
|
| 1327 |
|
To carry out isobaric-isothermal ensemble calculations {\sc oopse} |
| 1329 |
|
equations of motion,\cite{melchionna93} |
| 1330 |
|
|
| 1331 |
|
\begin{eqnarray} |
| 1332 |
< |
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right) \\ |
| 1333 |
< |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v} \\ |
| 1332 |
> |
\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1333 |
> |
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v}, \\ |
| 1334 |
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1335 |
< |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 1335 |
> |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right),\\ |
| 1336 |
|
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1337 |
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1338 |
< |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\ |
| 1338 |
> |
V}{\partial \mathsf{A}} \right) - \chi {\bf j}, \\ |
| 1339 |
|
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1340 |
< |
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\ |
| 1340 |
> |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1341 |
|
\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P - |
| 1342 |
< |
P_{\mathrm{target}} \right) \\ |
| 1343 |
< |
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta |
| 1342 |
> |
P_{\mathrm{target}} \right), \\ |
| 1343 |
> |
\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta . |
| 1344 |
|
\label{eq:melchionna1} |
| 1345 |
|
\end{eqnarray} |
| 1346 |
|
|
| 1353 |
|
volume can be calculated from the determinant of the matrix which |
| 1354 |
|
describes the box shape: |
| 1355 |
|
\begin{equation} |
| 1356 |
< |
\mathcal{V} = \det(\mathsf{H}) |
| 1356 |
> |
\mathcal{V} = \det(\mathsf{H}). |
| 1357 |
|
\end{equation} |
| 1358 |
|
|
| 1359 |
|
The NPTi integrator requires an instantaneous pressure. This quantity |
| 1361 |
|
\begin{equation} |
| 1362 |
|
\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{\mathcal{V}(t)} \left( |
| 1363 |
|
\sum_{i=1}^{N} m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) + |
| 1364 |
< |
\overleftrightarrow{\mathsf{W}}(t) |
| 1364 |
> |
\overleftrightarrow{\mathsf{W}}(t). |
| 1365 |
|
\end{equation} |
| 1366 |
|
The kinetic contribution to the pressure tensor utilizes the {\it |
| 1367 |
|
outer} product of the velocities denoted by the $\otimes$ symbol. The |
| 1370 |
|
r}_i$) with the forces between the same two atoms, |
| 1371 |
|
\begin{equation} |
| 1372 |
|
\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t) |
| 1373 |
< |
\otimes {\bf f}_{ij}(t) |
| 1373 |
> |
\otimes {\bf f}_{ij}(t). |
| 1374 |
|
\end{equation} |
| 1375 |
|
The instantaneous pressure is then simply obtained from the trace of |
| 1376 |
|
the Pressure tensor, |
| 1377 |
|
\begin{equation} |
| 1378 |
< |
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 1378 |
> |
P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t). |
| 1379 |
|
\right) |
| 1380 |
|
\end{equation} |
| 1381 |
|
|
| 1389 |
|
velocity-Verlet style 2 part algorithm: |
| 1390 |
|
|
| 1391 |
|
{\tt moveA:} |
| 1392 |
< |
\begin{eqnarray} |
| 1393 |
< |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 1394 |
< |
P(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\ |
| 1395 |
< |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 1396 |
< |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1397 |
< |
\left(\chi(t) + \eta(t) \right) \right) \\ |
| 1398 |
< |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 1399 |
< |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1400 |
< |
\chi(t) \right) \\ |
| 1401 |
< |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 1402 |
< |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\ |
| 1403 |
< |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 1404 |
< |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1405 |
< |
\right) \\ |
| 1406 |
< |
\eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B |
| 1407 |
< |
T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\ |
| 1408 |
< |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf |
| 1409 |
< |
v}\left(t + \delta t / 2 \right) + \eta(t + \delta t / 2)\left[ {\bf |
| 1410 |
< |
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\ |
| 1411 |
< |
\mathsf{H}(t + \delta t) & \leftarrow & e^{-\delta t \eta(t + \delta t |
| 1412 |
< |
/ 2)} \mathsf{H}(t) |
| 1413 |
< |
\end{eqnarray} |
| 1392 |
> |
\begin{align*} |
| 1393 |
> |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1394 |
> |
% |
| 1395 |
> |
P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} ,\\ |
| 1396 |
> |
% |
| 1397 |
> |
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1398 |
> |
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t) |
| 1399 |
> |
\left(\chi(t) + \eta(t) \right) \right), \\ |
| 1400 |
> |
% |
| 1401 |
> |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1402 |
> |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1403 |
> |
\chi(t) \right), \\ |
| 1404 |
> |
% |
| 1405 |
> |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1406 |
> |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1407 |
> |
\right) ,\\ |
| 1408 |
> |
% |
| 1409 |
> |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1410 |
> |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1411 |
> |
\right) ,\\ |
| 1412 |
> |
% |
| 1413 |
> |
\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h |
| 1414 |
> |
\mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t) |
| 1415 |
> |
- P_{\mathrm{target}} \right), \\ |
| 1416 |
> |
% |
| 1417 |
> |
{\bf r}(t + h) &\leftarrow {\bf r}(t) + h |
| 1418 |
> |
\left\{ {\bf v}\left(t + h / 2 \right) |
| 1419 |
> |
+ \eta(t + h / 2)\left[ {\bf r}(t + h) |
| 1420 |
> |
- {\bf R}_0 \right] \right\} ,\\ |
| 1421 |
> |
% |
| 1422 |
> |
\mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)} |
| 1423 |
> |
\mathsf{H}(t). |
| 1424 |
> |
\end{align*} |
| 1425 |
|
|
| 1426 |
|
Most of these equations are identical to their counterparts in the NVT |
| 1427 |
< |
integrator, but the propagation of positions to time $t + \delta t$ |
| 1427 |
> |
integrator, but the propagation of positions to time $t + h$ |
| 1428 |
|
depends on the positions at the same time. {\sc oopse} carries out |
| 1429 |
|
this step iteratively (with a limit of 5 passes through the iterative |
| 1430 |
|
loop). Also, the simulation box $\mathsf{H}$ is scaled uniformly for |
| 1431 |
|
one full time step by an exponential factor that depends on the value |
| 1432 |
|
of $\eta$ at time $t + |
| 1433 |
< |
\delta t / 2$. Reshaping the box uniformly also scales the volume of |
| 1433 |
> |
h / 2$. Reshaping the box uniformly also scales the volume of |
| 1434 |
|
the box by |
| 1435 |
|
\begin{equation} |
| 1436 |
< |
\mathcal{V}(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)} |
| 1436 |
> |
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}. |
| 1437 |
|
\mathcal{V}(t) |
| 1438 |
|
\end{equation} |
| 1439 |
|
|
| 1443 |
|
the same time value. |
| 1444 |
|
|
| 1445 |
|
{\tt moveB:} |
| 1446 |
< |
\begin{eqnarray} |
| 1447 |
< |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 1448 |
< |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 1449 |
< |
P(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\}, |
| 1450 |
< |
\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\ |
| 1451 |
< |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 1452 |
< |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 1453 |
< |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 1454 |
< |
\eta(t + \delta t) & \leftarrow & \eta(t + \delta t / 2) + |
| 1455 |
< |
\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2} |
| 1456 |
< |
\left( P(t + \delta t) - P_{\mathrm{target}} |
| 1457 |
< |
\right) \\ |
| 1458 |
< |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1459 |
< |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 1460 |
< |
\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t) |
| 1461 |
< |
(\chi(t + \delta t) + \eta(t + \delta t)) \right) \\ |
| 1462 |
< |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1463 |
< |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 1464 |
< |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 1465 |
< |
\chi(t + \delta t) \right) |
| 1466 |
< |
\end{eqnarray} |
| 1446 |
> |
\begin{align*} |
| 1447 |
> |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1448 |
> |
\left\{{\bf j}(t + h)\right\} ,\\ |
| 1449 |
> |
% |
| 1450 |
> |
P(t + h) &\leftarrow \left\{{\bf r}(t + h)\right\}, |
| 1451 |
> |
\left\{{\bf v}(t + h)\right\}, \\ |
| 1452 |
> |
% |
| 1453 |
> |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1454 |
> |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)} |
| 1455 |
> |
{T_{\mathrm{target}}} - 1 \right), \\ |
| 1456 |
> |
% |
| 1457 |
> |
\eta(t + h) &\leftarrow \eta(t + h / 2) + |
| 1458 |
> |
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1459 |
> |
\tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right), \\ |
| 1460 |
> |
% |
| 1461 |
> |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t |
| 1462 |
> |
+ h / 2 \right) + \frac{h}{2} \left( |
| 1463 |
> |
\frac{{\bf f}(t + h)}{m} - {\bf v}(t + h) |
| 1464 |
> |
(\chi(t + h) + \eta(t + h)) \right) ,\\ |
| 1465 |
> |
% |
| 1466 |
> |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
| 1467 |
> |
+ h / 2 \right) + \frac{h}{2} \left( {\bf |
| 1468 |
> |
\tau}^b(t + h) - {\bf j}(t + h) |
| 1469 |
> |
\chi(t + h) \right) . |
| 1470 |
> |
\end{align*} |
| 1471 |
|
|
| 1472 |
< |
Once again, since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ |
| 1473 |
< |
are required to caclculate $T(t + \delta t)$, $P(t + \delta t)$, $\chi(t + |
| 1474 |
< |
\delta t)$, and $\eta(t + \delta t)$, they indirectly depend on their |
| 1475 |
< |
own values at time $t + \delta t$. {\tt moveB} is therefore done in |
| 1476 |
< |
an iterative fashion until $\chi(t + \delta t)$ and $\eta(t + \delta |
| 1477 |
< |
t)$ become self-consistent. The relative tolerance for the |
| 1478 |
< |
self-consistency check defaults to a value of $\mbox{10}^{-6}$, but |
| 1457 |
< |
{\sc oopse} will terminate the iteration after 4 loops even if the |
| 1472 |
> |
Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required |
| 1473 |
> |
to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t + |
| 1474 |
> |
h)$, they indirectly depend on their own values at time $t + h$. {\tt |
| 1475 |
> |
moveB} is therefore done in an iterative fashion until $\chi(t + h)$ |
| 1476 |
> |
and $\eta(t + h)$ become self-consistent. The relative tolerance for |
| 1477 |
> |
the self-consistency check defaults to a value of $\mbox{10}^{-6}$, |
| 1478 |
> |
but {\sc oopse} will terminate the iteration after 4 loops even if the |
| 1479 |
|
consistency check has not been satisfied. |
| 1480 |
|
|
| 1481 |
|
The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is |
| 1502 |
|
{\tt moveB} portions of the algorithm. Details on the constraint |
| 1503 |
|
algorithms are given in section \ref{oopseSec:rattle}. |
| 1504 |
|
|
| 1505 |
< |
\subsubsection{\label{sec:NPTf}Constant-pressure integration with a |
| 1505 |
> |
\subsection{\label{sec:NPTf}Constant-pressure integration with a |
| 1506 |
|
flexible box (NPTf)} |
| 1507 |
|
|
| 1508 |
|
There is a relatively simple generalization of the |
| 1512 |
|
extended variables ($\overleftrightarrow{\eta}$) to control changes to |
| 1513 |
|
the box shape. The equations of motion for this method are |
| 1514 |
|
\begin{eqnarray} |
| 1515 |
< |
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right) \\ |
| 1515 |
> |
\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right), \\ |
| 1516 |
|
\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} + |
| 1517 |
< |
\chi \mathsf{1}) {\bf v} \\ |
| 1517 |
> |
\chi \cdot \mathsf{1}) {\bf v}, \\ |
| 1518 |
|
\dot{\mathsf{A}} & = & \mathsf{A} \cdot |
| 1519 |
< |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\ |
| 1519 |
> |
\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) ,\\ |
| 1520 |
|
\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1} |
| 1521 |
|
\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial |
| 1522 |
< |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\ |
| 1522 |
> |
V}{\partial \mathsf{A}} \right) - \chi {\bf j} ,\\ |
| 1523 |
|
\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left( |
| 1524 |
< |
\frac{T}{T_{\mathrm{target}}} - 1 \right) \\ |
| 1525 |
< |
\dot{\overleftrightarrow{eta}} & = & \frac{1}{\tau_{B}^2 f k_B |
| 1526 |
< |
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) \\ |
| 1527 |
< |
\dot{\mathsf{H}} & = & \overleftrightarrow{\eta} \cdot \mathsf{H} |
| 1524 |
> |
\frac{T}{T_{\mathrm{target}}} - 1 \right) ,\\ |
| 1525 |
> |
\dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B |
| 1526 |
> |
T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) ,\\ |
| 1527 |
> |
\dot{\mathsf{H}} & = & \overleftrightarrow{\eta} \cdot \mathsf{H} . |
| 1528 |
|
\label{eq:melchionna2} |
| 1529 |
|
\end{eqnarray} |
| 1530 |
|
|
| 1536 |
|
NPTi integration: |
| 1537 |
|
|
| 1538 |
|
{\tt moveA:} |
| 1539 |
< |
\begin{eqnarray} |
| 1540 |
< |
T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 1541 |
< |
\overleftrightarrow{\mathsf{P}}(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\ |
| 1542 |
< |
{\bf v}\left(t + \delta t / 2\right) & \leftarrow & {\bf |
| 1543 |
< |
v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - |
| 1544 |
< |
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
| 1545 |
< |
{\bf v}(t) \right) \\ |
| 1546 |
< |
{\bf j}\left(t + \delta t / 2 \right) & \leftarrow & {\bf |
| 1547 |
< |
j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1548 |
< |
\chi(t) \right) \\ |
| 1549 |
< |
\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t * |
| 1550 |
< |
{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\ |
| 1551 |
< |
\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) + |
| 1552 |
< |
\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1 |
| 1553 |
< |
\right) \\ |
| 1554 |
< |
\overleftrightarrow{\eta}(t + \delta t / 2) & \leftarrow & \overleftrightarrow{\eta}(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B |
| 1555 |
< |
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) - P_{\mathrm{target}}\mathsf{1} \right) \\ |
| 1556 |
< |
{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf |
| 1557 |
< |
v}\left(t + \delta t / 2 \right) + \overleftrightarrow{\eta}(t + |
| 1558 |
< |
\delta t / 2) \cdot \left[ {\bf |
| 1559 |
< |
r}(t + \delta t) - {\bf R}_0 \right] \right\} \\ |
| 1560 |
< |
\mathsf{H}(t + \delta t) & \leftarrow & \mathsf{H}(t) \cdot e^{-\delta t |
| 1561 |
< |
\overleftrightarrow{\eta}(t + \delta t / 2)} |
| 1562 |
< |
\end{eqnarray} |
| 1539 |
> |
\begin{align*} |
| 1540 |
> |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\ |
| 1541 |
> |
% |
| 1542 |
> |
\overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\}, |
| 1543 |
> |
\left\{{\bf v}(t)\right\} ,\\ |
| 1544 |
> |
% |
| 1545 |
> |
{\bf v}\left(t + h / 2\right) &\leftarrow {\bf v}(t) |
| 1546 |
> |
+ \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - |
| 1547 |
> |
\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot |
| 1548 |
> |
{\bf v}(t) \right), \\ |
| 1549 |
> |
% |
| 1550 |
> |
{\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) |
| 1551 |
> |
+ \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t) |
| 1552 |
> |
\chi(t) \right), \\ |
| 1553 |
> |
% |
| 1554 |
> |
\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h * |
| 1555 |
> |
{\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1} |
| 1556 |
> |
\right), \\ |
| 1557 |
> |
% |
| 1558 |
> |
\chi\left(t + h / 2 \right) &\leftarrow \chi(t) + |
| 1559 |
> |
\frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} |
| 1560 |
> |
- 1 \right), \\ |
| 1561 |
> |
% |
| 1562 |
> |
\overleftrightarrow{\eta}(t + h / 2) &\leftarrow |
| 1563 |
> |
\overleftrightarrow{\eta}(t) + \frac{h \mathcal{V}(t)}{2 N k_B |
| 1564 |
> |
T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) |
| 1565 |
> |
- P_{\mathrm{target}}\mathsf{1} \right), \\ |
| 1566 |
> |
% |
| 1567 |
> |
{\bf r}(t + h) &\leftarrow {\bf r}(t) + h \left\{ {\bf v} |
| 1568 |
> |
\left(t + h / 2 \right) + \overleftrightarrow{\eta}(t + |
| 1569 |
> |
h / 2) \cdot \left[ {\bf r}(t + h) |
| 1570 |
> |
- {\bf R}_0 \right] \right\}, \\ |
| 1571 |
> |
% |
| 1572 |
> |
\mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h |
| 1573 |
> |
\overleftrightarrow{\eta}(t + h / 2)} . |
| 1574 |
> |
\end{align*} |
| 1575 |
|
{\sc oopse} uses a power series expansion truncated at second order |
| 1576 |
|
for the exponential operation which scales the simulation box. |
| 1577 |
|
|
| 1579 |
|
NPTi integrator: |
| 1580 |
|
|
| 1581 |
|
{\tt moveB:} |
| 1582 |
< |
\begin{eqnarray} |
| 1583 |
< |
T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\}, |
| 1584 |
< |
\left\{{\bf j}(t + \delta t)\right\} \\ |
| 1585 |
< |
\overleftrightarrow{\mathsf{P}}(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\}, |
| 1586 |
< |
\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\ |
| 1587 |
< |
\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t / |
| 1588 |
< |
2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta |
| 1589 |
< |
t)}{T_{\mathrm{target}}} - 1 \right) \\ |
| 1590 |
< |
\overleftrightarrow{\eta}(t + \delta t) & \leftarrow & \overleftrightarrow{\eta}(t + \delta t / 2) + |
| 1591 |
< |
\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2} |
| 1592 |
< |
\left( \overleftrightarrow{P}(t + \delta t) - P_{\mathrm{target}}\mathsf{1} |
| 1593 |
< |
\right) \\ |
| 1594 |
< |
{\bf v}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1595 |
< |
v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( |
| 1596 |
< |
\frac{{\bf f}(t + \delta t)}{m} - |
| 1597 |
< |
(\chi(t + \delta t)\mathsf{1} + \overleftrightarrow{\eta}(t + \delta |
| 1598 |
< |
t)) \right) \cdot {\bf v}(t + \delta t) \\ |
| 1599 |
< |
{\bf j}\left(t + \delta t \right) & \leftarrow & {\bf |
| 1600 |
< |
j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf |
| 1601 |
< |
\tau}^b(t + \delta t) - {\bf j}(t + \delta t) |
| 1602 |
< |
\chi(t + \delta t) \right) |
| 1603 |
< |
\end{eqnarray} |
| 1582 |
> |
\begin{align*} |
| 1583 |
> |
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1584 |
> |
\left\{{\bf j}(t + h)\right\}, \\ |
| 1585 |
> |
% |
| 1586 |
> |
\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r} |
| 1587 |
> |
(t + h)\right\}, \left\{{\bf v}(t |
| 1588 |
> |
+ h)\right\}, \left\{{\bf f}(t + h)\right\} ,\\ |
| 1589 |
> |
% |
| 1590 |
> |
\chi\left(t + h \right) &\leftarrow \chi\left(t + h / |
| 1591 |
> |
2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+ |
| 1592 |
> |
h)}{T_{\mathrm{target}}} - 1 \right), \\ |
| 1593 |
> |
% |
| 1594 |
> |
\overleftrightarrow{\eta}(t + h) &\leftarrow |
| 1595 |
> |
\overleftrightarrow{\eta}(t + h / 2) + |
| 1596 |
> |
\frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h) |
| 1597 |
> |
\tau_B^2} \left( \overleftrightarrow{P}(t + h) |
| 1598 |
> |
- P_{\mathrm{target}}\mathsf{1} \right) ,\\ |
| 1599 |
> |
% |
| 1600 |
> |
{\bf v}\left(t + h \right) &\leftarrow {\bf v}\left(t |
| 1601 |
> |
+ h / 2 \right) + \frac{h}{2} \left( |
| 1602 |
> |
\frac{{\bf f}(t + h)}{m} - |
| 1603 |
> |
(\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t |
| 1604 |
> |
+ h)) \right) \cdot {\bf v}(t + h), \\ |
| 1605 |
> |
% |
| 1606 |
> |
{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t |
| 1607 |
> |
+ h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t |
| 1608 |
> |
+ h) - {\bf j}(t + h) \chi(t + h) \right) . |
| 1609 |
> |
\end{align*} |
| 1610 |
|
|
| 1611 |
|
The iterative schemes for both {\tt moveA} and {\tt moveB} are |
| 1612 |
|
identical to those described for the NPTi integrator. |
| 1624 |
|
simulations. Liquids have very small restoring forces in the |
| 1625 |
|
off-diagonal directions, and the simulation box can very quickly form |
| 1626 |
|
elongated and sheared geometries which become smaller than the |
| 1627 |
< |
electrostatic or Lennard-Jones cutoff radii. It finds most use in |
| 1628 |
< |
simulating crystals or liquid crystals which assume non-orthorhombic |
| 1629 |
< |
geometries. |
| 1627 |
> |
electrostatic or Lennard-Jones cutoff radii. The NPTf integrator |
| 1628 |
> |
finds most use in simulating crystals or liquid crystals which assume |
| 1629 |
> |
non-orthorhombic geometries. |
| 1630 |
|
|
| 1631 |
< |
\subsubsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)} |
| 1631 |
> |
\subsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)} |
| 1632 |
|
|
| 1633 |
|
There is one additional extended system integrator which is somewhat |
| 1634 |
|
simpler than the NPTf method described above. In this case, the three |
| 1656 |
|
oopse}, we have implemented the {\sc rattle} algorithm of |
| 1657 |
|
Andersen.\cite{andersen83} The algorithm is a velocity verlet |
| 1658 |
|
formulation of the {\sc shake} method\cite{ryckaert77} of iteratively |
| 1659 |
< |
solving the Lagrange multipliers of constraint. The system of lagrange |
| 1659 |
> |
solving the Lagrange multipliers of constraint. The system of Lagrange |
| 1660 |
|
multipliers allows one to reformulate the equations of motion with |
| 1661 |
|
explicit constraint forces.\cite{fowles99:lagrange} |
| 1662 |
|
|
| 1671 |
|
\delta\int_{t_1}^{t_2}L\, dt = |
| 1672 |
|
\int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i} |
| 1673 |
|
- \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i} |
| 1674 |
< |
\biggr ) \biggr] \delta q_i \, dt = 0 |
| 1674 |
> |
\biggr ) \biggr] \delta q_i \, dt = 0. |
| 1675 |
|
\label{oopseEq:lm2} |
| 1676 |
|
\end{equation} |
| 1677 |
|
Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$ |
| 1679 |
|
instant of time, giving: |
| 1680 |
|
\begin{align} |
| 1681 |
|
\delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 % |
| 1682 |
< |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 \\ |
| 1682 |
> |
+ \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 ,\\ |
| 1683 |
|
% |
| 1684 |
|
\frac{\partial\sigma}{\partial q_1} \delta q_1 &= % |
| 1685 |
< |
- \frac{\partial\sigma}{\partial q_2} \delta q_2 \\ |
| 1685 |
> |
- \frac{\partial\sigma}{\partial q_2} \delta q_2, \\ |
| 1686 |
|
% |
| 1687 |
|
\delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1688 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1 |
| 1688 |
> |
\frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1. |
| 1689 |
|
\end{align} |
| 1690 |
|
Substituted back into Eq.~\ref{oopseEq:lm2}, |
| 1691 |
|
\begin{equation} |
| 1695 |
|
- \biggl( \frac{\partial L}{\partial q_1} |
| 1696 |
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1697 |
|
\biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / % |
| 1698 |
< |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0 |
| 1698 |
> |
\frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0. |
| 1699 |
|
\label{oopseEq:lm3} |
| 1700 |
|
\end{equation} |
| 1701 |
|
Leading to, |
| 1705 |
|
\biggr)}{\frac{\partial\sigma}{\partial q_1}} = |
| 1706 |
|
\frac{\biggl(\frac{\partial L}{\partial q_2} |
| 1707 |
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2} |
| 1708 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_2}} |
| 1708 |
> |
\biggr)}{\frac{\partial\sigma}{\partial q_2}}. |
| 1709 |
|
\label{oopseEq:lm4} |
| 1710 |
|
\end{equation} |
| 1711 |
|
This relation can only be statisfied, if both are equal to a single |
| 1713 |
|
\begin{align} |
| 1714 |
|
\frac{\biggl(\frac{\partial L}{\partial q_1} |
| 1715 |
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1716 |
< |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t) \\ |
| 1716 |
> |
\biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t), \\ |
| 1717 |
|
% |
| 1718 |
|
\frac{\partial L}{\partial q_1} |
| 1719 |
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &= |
| 1720 |
< |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} \\ |
| 1720 |
> |
-\lambda(t)\,\frac{\partial\sigma}{\partial q_1} ,\\ |
| 1721 |
|
% |
| 1722 |
|
\frac{\partial L}{\partial q_1} |
| 1723 |
|
- \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} |
| 1724 |
< |
+ \mathcal{G}_i &= 0 |
| 1724 |
> |
+ \mathcal{G}_i &= 0, |
| 1725 |
|
\end{align} |
| 1726 |
< |
Where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1726 |
> |
where $\mathcal{G}_i$, the force of constraint on $i$, is: |
| 1727 |
|
\begin{equation} |
| 1728 |
< |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1} |
| 1728 |
> |
\mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}. |
| 1729 |
|
\label{oopseEq:lm5} |
| 1730 |
|
\end{equation} |
| 1731 |
|
|
| 1751 |
|
\begin{align} |
| 1752 |
|
\sigma_{ij}[\mathbf{r}(t)] \equiv |
| 1753 |
|
[ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2 - d_{ij}^2 &= 0 % |
| 1754 |
< |
\label{oopseEq:c1} \\ |
| 1754 |
> |
\label{oopseEq:c1}, \\ |
| 1755 |
|
% |
| 1756 |
|
[\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot |
| 1757 |
< |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 \label{oopseEq:c2} |
| 1757 |
> |
[\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 .\label{oopseEq:c2} |
| 1758 |
|
\end{align} |
| 1759 |
|
Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is |
| 1760 |
|
the constrained distance between atom $i$ and |
| 1762 |
|
be perpendicular to the bond vector, so that the bond can neither grow |
| 1763 |
|
nor shrink. The constrained dynamics equations become: |
| 1764 |
|
\begin{equation} |
| 1765 |
< |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i |
| 1765 |
> |
m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i, |
| 1766 |
|
\label{oopseEq:r1} |
| 1767 |
|
\end{equation} |
| 1768 |
< |
Where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$, |
| 1768 |
> |
where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$, |
| 1769 |
|
and are defined: |
| 1770 |
|
\begin{equation} |
| 1771 |
< |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij} |
| 1771 |
> |
\mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}. |
| 1772 |
|
\label{oopseEq:r2} |
| 1773 |
|
\end{equation} |
| 1774 |
|
|
| 1777 |
|
\mathbf{r}_i(t+h) &= |
| 1778 |
|
\mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) + |
| 1779 |
|
\frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) + |
| 1780 |
< |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1} \\ |
| 1780 |
> |
\mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1}, \\ |
| 1781 |
|
% |
| 1782 |
|
\mathbf{\dot{r}}_i(t+h) &= |
| 1783 |
|
\mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} |
| 1784 |
|
\Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) + |
| 1785 |
< |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] % |
| 1785 |
> |
\mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] ,% |
| 1786 |
|
\label{oopseEq:vv2} |
| 1787 |
|
\end{align} |
| 1788 |
< |
Where: |
| 1788 |
> |
where: |
| 1789 |
|
\begin{align} |
| 1790 |
|
\mathbf{\mathcal{G}}_{Ri}(t) &= |
| 1791 |
< |
-2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) \\ |
| 1791 |
> |
-2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) ,\\ |
| 1792 |
|
% |
| 1793 |
|
\mathbf{\mathcal{G}}_{Vi}(t+h) &= |
| 1794 |
< |
-2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h) |
| 1794 |
> |
-2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h). |
| 1795 |
|
\end{align} |
| 1796 |
|
Next, define: |
| 1797 |
|
\begin{align} |
| 1798 |
< |
g_{ij} &= h \lambda_{Rij}(t) \\ |
| 1799 |
< |
k_{ij} &= h \lambda_{Vij}(t+h) \\ |
| 1798 |
> |
g_{ij} &= h \lambda_{Rij}(t) ,\\ |
| 1799 |
> |
k_{ij} &= h \lambda_{Vij}(t+h), \\ |
| 1800 |
|
\mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t) |
| 1801 |
< |
- \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t) |
| 1801 |
> |
- \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t). |
| 1802 |
|
\end{align} |
| 1803 |
|
Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2} |
| 1804 |
|
can be rewritten as, |
| 1805 |
|
\begin{align} |
| 1806 |
< |
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i \\ |
| 1806 |
> |
\mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i ,\\ |
| 1807 |
|
% |
| 1808 |
|
\mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1809 |
< |
-\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h) |
| 1809 |
> |
-\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h). |
| 1810 |
|
\end{align} |
| 1811 |
|
|
| 1812 |
|
To integrate the equations of motion, the {\sc rattle} algorithm first |
| 1813 |
|
solves for $\mathbf{r}(t+h)$. Let, |
| 1814 |
|
\begin{equation} |
| 1815 |
< |
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t) |
| 1815 |
> |
\mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t). |
| 1816 |
|
\end{equation} |
| 1817 |
|
Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next |
| 1818 |
|
pick a constraint $j$, and let, |
| 1819 |
|
\begin{equation} |
| 1820 |
|
\mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t) |
| 1821 |
< |
- \mathbf{r}_j(t) + h\mathbf{q}_j(t) |
| 1821 |
> |
- \mathbf{r}_j(t) + h\mathbf{q}_j(t). |
| 1822 |
|
\label{oopseEq:ra1} |
| 1823 |
|
\end{equation} |
| 1824 |
|
If |
| 1829 |
|
positions. First we define a test corrected configuration as, |
| 1830 |
|
\begin{align} |
| 1831 |
|
\mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i - |
| 1832 |
< |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] \\ |
| 1832 |
> |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] ,\\ |
| 1833 |
|
% |
| 1834 |
|
\mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j + |
| 1835 |
< |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr] |
| 1835 |
> |
g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr]. |
| 1836 |
|
\end{align} |
| 1837 |
|
And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2 |
| 1838 |
|
= d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the |
| 1839 |
|
approximation, |
| 1840 |
|
\begin{equation} |
| 1841 |
|
g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)] |
| 1842 |
< |
(\frac{1}{m_i} + \frac{1}{m_j})} |
| 1842 |
> |
(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1843 |
|
\end{equation} |
| 1844 |
|
Although not an exact solution for $g_{ij}$, as this is an iterative |
| 1845 |
|
scheme overall, the eventual solution will converge. With a trial |
| 1846 |
|
$g_{ij}$, the new $\mathbf{q}$'s become, |
| 1847 |
|
\begin{align} |
| 1848 |
|
\mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\, |
| 1849 |
< |
\frac{\mathbf{r}_{ij}(t)}{m_i} \\ |
| 1849 |
> |
\frac{\mathbf{r}_{ij}(t)}{m_i} ,\\ |
| 1850 |
|
% |
| 1851 |
|
\mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\, |
| 1852 |
< |
\frac{\mathbf{r}_{ij}(t)}{m_j} |
| 1852 |
> |
\frac{\mathbf{r}_{ij}(t)}{m_j} . |
| 1853 |
|
\end{align} |
| 1854 |
|
The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until |
| 1855 |
|
all constraints are satisfied. |
| 1857 |
|
The second step of {\sc rattle}, is to then update the velocities. The |
| 1858 |
|
step starts with, |
| 1859 |
|
\begin{equation} |
| 1860 |
< |
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h) |
| 1860 |
> |
\mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h). |
| 1861 |
|
\end{equation} |
| 1862 |
|
Next we pick a constraint $j$, and calculate the dot product $\ell$. |
| 1863 |
|
\begin{equation} |
| 1864 |
< |
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h) |
| 1864 |
> |
\ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h). |
| 1865 |
|
\label{oopseEq:rv1} |
| 1866 |
|
\end{equation} |
| 1867 |
|
Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be |
| 1869 |
|
corrections are made to the $i$ and $j$ velocities. |
| 1870 |
|
\begin{align} |
| 1871 |
|
\mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij} |
| 1872 |
< |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i} \\ |
| 1872 |
> |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i}, \\ |
| 1873 |
|
% |
| 1874 |
|
\mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij} |
| 1875 |
< |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j} |
| 1875 |
> |
\frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j}. |
| 1876 |
|
\end{align} |
| 1877 |
|
Like in the previous step, we select a value for $k_{ij}$ such that |
| 1878 |
|
$\ell$ is zero. |
| 1879 |
|
\begin{equation} |
| 1880 |
< |
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})} |
| 1880 |
> |
k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})}. |
| 1881 |
|
\end{equation} |
| 1882 |
|
The test velocities, $\mathbf{\dot{r}}^T_i$ and |
| 1883 |
|
$\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and |
| 1893 |
|
coefficient can be calculated from the deviation of the instantaneous |
| 1894 |
|
force from its mean force. |
| 1895 |
|
\begin{equation} |
| 1896 |
< |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T |
| 1896 |
> |
\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T, |
| 1897 |
|
\end{equation} |
| 1898 |
|
where% |
| 1899 |
|
\begin{equation} |
| 1900 |
< |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle |
| 1900 |
> |
\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle. |
| 1901 |
|
\end{equation} |
| 1902 |
|
|
| 1903 |
|
|
| 1904 |
|
If the time-dependent friction decays rapidly, the static friction |
| 1905 |
|
coefficient can be approximated by |
| 1906 |
|
\begin{equation} |
| 1907 |
< |
\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt |
| 1907 |
> |
\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt. |
| 1908 |
|
\end{equation} |
| 1909 |
|
Allowing diffusion constant to then be calculated through the |
| 1910 |
|
Einstein relation:\cite{Marrink94} |
| 1911 |
|
\begin{equation} |
| 1912 |
|
D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty |
| 1913 |
< |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}% |
| 1913 |
> |
}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.% |
| 1914 |
|
\end{equation} |
| 1915 |
|
|
| 1916 |
|
The Z-Constraint method, which fixes the z coordinates of the |
| 1925 |
|
|
| 1926 |
|
After the force calculation, define $G_\alpha$ as |
| 1927 |
|
\begin{equation} |
| 1928 |
< |
G_{\alpha} = \sum_i F_{\alpha i} |
| 1928 |
> |
G_{\alpha} = \sum_i F_{\alpha i}, |
| 1929 |
|
\label{oopseEq:zc1} |
| 1930 |
|
\end{equation} |
| 1931 |
< |
Where $F_{\alpha i}$ is the force in the z direction of atom $i$ in |
| 1931 |
> |
where $F_{\alpha i}$ is the force in the z direction of atom $i$ in |
| 1932 |
|
z-constrained molecule $\alpha$. The forces of the z constrained |
| 1933 |
|
molecule are then set to: |
| 1934 |
|
\begin{equation} |
| 1935 |
|
F_{\alpha i} = F_{\alpha i} - |
| 1936 |
< |
\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}} |
| 1936 |
> |
\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}. |
| 1937 |
|
\end{equation} |
| 1938 |
|
Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained |
| 1939 |
|
molecule. Having rescaled the forces, the velocities must also be |
| 1941 |
|
direction. |
| 1942 |
|
\begin{equation} |
| 1943 |
|
v_{\alpha i} = v_{\alpha i} - |
| 1944 |
< |
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}} |
| 1944 |
> |
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}, |
| 1945 |
|
\end{equation} |
| 1946 |
< |
Where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction. |
| 1946 |
> |
where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction. |
| 1947 |
|
Lastly, all of the accumulated z constrained forces must be subtracted |
| 1948 |
|
from the system to keep the system center of mass from drifting. |
| 1949 |
|
\begin{equation} |
| 1950 |
|
F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}} |
| 1951 |
< |
{\sum_{\beta}\sum_i m_{\beta i}} |
| 1951 |
> |
{\sum_{\beta}\sum_i m_{\beta i}}, |
| 1952 |
|
\end{equation} |
| 1953 |
< |
Where $\beta$ are all of the unconstrained molecules in the system. |
| 1953 |
> |
where $\beta$ are all of the unconstrained molecules in the |
| 1954 |
> |
system. Similarly, the velocities of the unconstrained molecules must |
| 1955 |
> |
also be scaled. |
| 1956 |
> |
\begin{equation} |
| 1957 |
> |
v_{\beta i} = v_{\beta i} + \sum_{\alpha} |
| 1958 |
> |
\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}. |
| 1959 |
> |
\end{equation} |
| 1960 |
|
|
| 1961 |
|
At the very beginning of the simulation, the molecules may not be at their |
| 1962 |
|
constrained positions. To move a z-constrained molecule to its specified |
| 1963 |
|
position, a simple harmonic potential is used |
| 1964 |
|
\begin{equation} |
| 1965 |
< |
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2}% |
| 1965 |
> |
U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},% |
| 1966 |
|
\end{equation} |
| 1967 |
|
where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$ is the |
| 1968 |
|
current $z$ coordinate of the center of mass of the constrained molecule, and |
| 1970 |
|
on the z-constrained molecule at time $t$ can be calculated by |
| 1971 |
|
\begin{equation} |
| 1972 |
|
F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}= |
| 1973 |
< |
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}) |
| 1973 |
> |
-k_{\text{Harmonic}}(z(t)-z_{\text{cons}}). |
| 1974 |
|
\end{equation} |
| 1975 |
|
|
| 1976 |
|
\section{\label{oopseSec:props}Trajectory Analysis} |
| 1984 |
|
can be found in Table~\ref{oopseTb:gofrs} |
| 1985 |
|
|
| 1986 |
|
\begin{table} |
| 1987 |
< |
\caption[The list of pair correlations in \texttt{staticProps}]{The different pair correlations in \texttt{staticProps} along with whether atom A or B must be directional.} |
| 1987 |
> |
\caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}} |
| 1988 |
|
\label{oopseTb:gofrs} |
| 1989 |
|
\begin{center} |
| 1990 |
|
\begin{tabular}{|l|c|c|} |
| 1997 |
|
$\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &% |
| 1998 |
|
both \\ \hline |
| 1999 |
|
\end{tabular} |
| 2000 |
+ |
\begin{minipage}{\linewidth} |
| 2001 |
+ |
\centering |
| 2002 |
+ |
\vspace{2mm} |
| 2003 |
+ |
The third column specifies which atom, if any, need be a directional entity. |
| 2004 |
+ |
\end{minipage} |
| 2005 |
|
\end{center} |
| 2006 |
|
\end{table} |
| 2007 |
|
|
| 2009 |
|
\begin{equation} |
| 2010 |
|
g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %% |
| 2011 |
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} %% |
| 2012 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle \label{eq:gofr} |
| 2012 |
> |
\delta( r - |\mathbf{r}_{ij}|) \rangle, \label{eq:gofr} |
| 2013 |
|
\end{equation} |
| 2014 |
< |
Where $\mathbf{r}_{ij}$ is the vector |
| 2014 |
> |
where $\mathbf{r}_{ij}$ is the vector |
| 2015 |
|
\begin{equation*} |
| 2016 |
< |
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i \notag |
| 2016 |
> |
\mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i, \notag |
| 2017 |
|
\end{equation*} |
| 2018 |
|
and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over |
| 2019 |
|
the expected pair density at a given $r$. |
| 2030 |
|
g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2031 |
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2032 |
|
\delta( \cos \theta - \cos \theta_{ij}) |
| 2033 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle |
| 2033 |
> |
\delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2034 |
|
\label{eq:gofrCosTheta} |
| 2035 |
|
\end{equation} |
| 2036 |
< |
Where |
| 2036 |
> |
Here |
| 2037 |
|
\begin{equation*} |
| 2038 |
< |
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij} |
| 2038 |
> |
\cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij}, |
| 2039 |
|
\end{equation*} |
| 2040 |
< |
Here $\mathbf{\hat{i}}$ is the unit directional vector of species $i$ |
| 2040 |
> |
where $\mathbf{\hat{i}}$ is the unit directional vector of species $i$ |
| 2041 |
|
and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector |
| 2042 |
|
$\mathbf{r}_{ij}$. |
| 2043 |
|
|
| 2047 |
|
\frac{V}{N_{\text{A}}N_{\text{B}}}\langle |
| 2048 |
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2049 |
|
\delta( \cos \omega - \cos \omega_{ij}) |
| 2050 |
< |
\delta( r - |\mathbf{r}_{ij}|) \rangle \label{eq:gofrCosOmega} |
| 2050 |
> |
\delta( r - |\mathbf{r}_{ij}|) \rangle. \label{eq:gofrCosOmega} |
| 2051 |
|
\end{equation} |
| 2052 |
|
Here |
| 2053 |
|
\begin{equation*} |
| 2054 |
< |
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}} |
| 2054 |
> |
\cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}}. |
| 2055 |
|
\end{equation*} |
| 2056 |
|
Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit |
| 2057 |
|
directional vectors of species $i$ and $j$. |
| 2064 |
|
\sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2065 |
|
\delta( x - x_{ij}) |
| 2066 |
|
\delta( y - y_{ij}) |
| 2067 |
< |
\delta( z - z_{ij}) \rangle |
| 2067 |
> |
\delta( z - z_{ij}) \rangle, |
| 2068 |
|
\end{equation} |
| 2069 |
< |
Where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$ |
| 2069 |
> |
where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$ |
| 2070 |
|
components respectively of vector $\mathbf{r}_{ij}$. |
| 2071 |
|
|
| 2072 |
|
The final pair correlation is similar to |
| 2075 |
|
\begin{equation}\label{eq:cosOmegaOfR} |
| 2076 |
|
\langle \cos \omega \rangle_{\text{AB}}(r) = |
| 2077 |
|
\langle \sum_{i \in \text{A}} \sum_{j \in \text{B}} |
| 2078 |
< |
(\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle |
| 2078 |
> |
(\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle. |
| 2079 |
|
\end{equation} |
| 2080 |
|
Here $\cos \omega_{ij}$ is defined in the same way as in |
| 2081 |
|
Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair |
| 2087 |
|
The dynamic properties of a trajectory are calculated with the program |
| 2088 |
|
\texttt{dynamicProps}. The program calculates the following properties: |
| 2089 |
|
\begin{gather} |
| 2090 |
< |
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle \label{eq:rms}\\ |
| 2091 |
< |
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle \label{eq:velCorr} \\ |
| 2092 |
< |
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle \label{eq:angularVelCorr} |
| 2090 |
> |
\langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle, \label{eq:rms}\\ |
| 2091 |
> |
\langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle, \label{eq:velCorr} \\ |
| 2092 |
> |
\langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle. \label{eq:angularVelCorr} |
| 2093 |
|
\end{gather} |
| 2094 |
|
|
| 2095 |
|
Eq.~\ref{eq:rms} is the root mean square displacement function. Which |
| 2098 |
|
coefficients because of the Einstein Relation, which is valid at long |
| 2099 |
|
times.\cite{allen87:csl} |
| 2100 |
|
\begin{equation} |
| 2101 |
< |
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle |
| 2101 |
> |
2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle. |
| 2102 |
|
\label{oopseEq:einstein} |
| 2103 |
|
\end{equation} |
| 2104 |
|
|
| 2170 |
|
Algorithmically simplest of the three methods is atomic decomposition |
| 2171 |
|
where N particles in a simulation are split among P processors for the |
| 2172 |
|
duration of the simulation. Computational cost scales as an optimal |
| 2173 |
< |
$O(N/P)$ for atomic decomposition. Unfortunately all processors must |
| 2174 |
< |
communicate positions and forces with all other processors at every |
| 2175 |
< |
force evaluation, leading communication costs to scale as an |
| 2176 |
< |
unfavorable $O(N)$, \emph{independent of the number of processors}. This |
| 2177 |
< |
communication bottleneck led to the development of spatial and force |
| 2178 |
< |
decomposition methods in which communication among processors scales |
| 2179 |
< |
much more favorably. Spatial or domain decomposition divides the |
| 2180 |
< |
physical spatial domain into 3D boxes in which each processor is |
| 2181 |
< |
responsible for calculation of forces and positions of particles |
| 2182 |
< |
located in its box. Particles are reassigned to different processors |
| 2183 |
< |
as they move through simulation space. To calculate forces on a given |
| 2184 |
< |
particle, a processor must know the positions of particles within some |
| 2185 |
< |
cutoff radius located on nearby processors instead of the positions of |
| 2186 |
< |
particles on all processors. Both communication between processors and |
| 2187 |
< |
computation scale as $O(N/P)$ in the spatial method. However, spatial |
| 2173 |
> |
$\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately all |
| 2174 |
> |
processors must communicate positions and forces with all other |
| 2175 |
> |
processors at every force evaluation, leading communication costs to |
| 2176 |
> |
scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the |
| 2177 |
> |
number of processors}. This communication bottleneck led to the |
| 2178 |
> |
development of spatial and force decomposition methods in which |
| 2179 |
> |
communication among processors scales much more favorably. Spatial or |
| 2180 |
> |
domain decomposition divides the physical spatial domain into 3D boxes |
| 2181 |
> |
in which each processor is responsible for calculation of forces and |
| 2182 |
> |
positions of particles located in its box. Particles are reassigned to |
| 2183 |
> |
different processors as they move through simulation space. To |
| 2184 |
> |
calculate forces on a given particle, a processor must know the |
| 2185 |
> |
positions of particles within some cutoff radius located on nearby |
| 2186 |
> |
processors instead of the positions of particles on all |
| 2187 |
> |
processors. Both communication between processors and computation |
| 2188 |
> |
scale as $\mathcal{O}(N/P)$ in the spatial method. However, spatial |
| 2189 |
|
decomposition adds algorithmic complexity to the simulation code and |
| 2190 |
|
is not very efficient for small N since the overall communication |
| 2191 |
< |
scales as the surface to volume ratio $(N/P)^{2/3}$ in three |
| 2192 |
< |
dimensions. |
| 2191 |
> |
scales as the surface to volume ratio $\mathcal{O}(N/P)^{2/3}$ in |
| 2192 |
> |
three dimensions. |
| 2193 |
|
|
| 2194 |
|
The parallelization method used in {\sc oopse} is the force |
| 2195 |
|
decomposition method. Force decomposition assigns particles to |
| 2198 |
|
and column processor groups. Forces are calculated on particles in a |
| 2199 |
|
given row by particles located in that processors column |
| 2200 |
|
assignment. Force decomposition is less complex to implement than the |
| 2201 |
< |
spatial method but still scales computationally as $O(N/P)$ and scales |
| 2202 |
< |
as $O(N/\sqrt{P})$ in communication cost. Plimpton has also found that |
| 2203 |
< |
force decompositions scale more favorably than spatial decompositions |
| 2204 |
< |
for systems up to 10,000 atoms and favorably compete with spatial |
| 2205 |
< |
methods up to 100,000 atoms.\cite{plimpton95} |
| 2201 |
> |
spatial method but still scales computationally as $\mathcal{O}(N/P)$ |
| 2202 |
> |
and scales as $\mathcal{O}(N/\sqrt{P})$ in communication |
| 2203 |
> |
cost. Plimpton has also found that force decompositions scale more |
| 2204 |
> |
favorably than spatial decompositions for systems up to 10,000 atoms |
| 2205 |
> |
and favorably compete with spatial methods up to 100,000 |
| 2206 |
> |
atoms.\cite{plimpton95} |
| 2207 |
|
|
| 2208 |
|
\subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis} |
| 2209 |
|
|
| 2229 |
|
\texttt{dynamicProps} will calculate all of the time correlation frame |
| 2230 |
|
pairs within the block. After in-block correlations are complete, a |
| 2231 |
|
second block of the trajectory is read, and the cross correlations are |
| 2232 |
< |
calculated between the two blocks. this second block is then freed and |
| 2232 |
> |
calculated between the two blocks. This second block is then freed and |
| 2233 |
|
then incremented and the process repeated until the end of the |
| 2234 |
|
trajectory. Once the end is reached, the first block is freed then |
| 2235 |
|
incremented, and the again the internal time correlations are |
| 2259 |
|
z-constraint method. |
| 2260 |
|
|
| 2261 |
|
These features are all brought together in a single open-source |
| 2262 |
< |
program. Allowing researchers to not only benefit from |
| 2262 |
> |
program. This allows researchers to not only benefit from |
| 2263 |
|
{\sc oopse}, but also contribute to {\sc oopse}'s development as |
| 2264 |
< |
well.Documentation and source code for {\sc oopse} can be downloaded |
| 2213 |
< |
from \texttt{http://www.openscience.org/oopse/}. |
| 2264 |
> |
well. |
| 2265 |
|
|
