| 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 |
| 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 an example \texttt{.bass} file that |
| 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 |
|
|
| 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 |
| 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 |
| 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]A portion of a \texttt{.bass} file showing a simulation including {\sc ssd} water.},label={sch:ssd}] |
| 581 |
|
|
| 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 |
| 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}_{\alpha}$ is the column vector of the $\alpha$ 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$ |
| 678 |
|
first convert it to its corresponding vector in box space, and then, |
| 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 |
| 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 |
| 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 |
| 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 |
|
$\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular |
| 914 |
|
momentum and moment of inertia tensor respectively, and the |
| 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, |
| 1003 |
|
{\tt moveA:} |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) |
| 1015 |
> |
(t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right). |
| 1016 |
|
\end{align*} |
| 1017 |
|
|
| 1018 |
|
In this context, the $\mathrm{rotate}$ function is the reversible product |
| 1020 |
|
\begin{equation} |
| 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} |
| 1057 |
|
{\tt doForces:} |
| 1058 |
|
\begin{align*} |
| 1059 |
|
{\bf f}(t + h) &\leftarrow |
| 1060 |
< |
- \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\ |
| 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}} \\ |
| 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) |
| 1066 |
> |
\cdot {\bf \tau}^s(t + h). |
| 1067 |
|
\end{align*} |
| 1068 |
|
|
| 1069 |
|
{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix |
| 1074 |
|
{\tt moveB:} |
| 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) \\ |
| 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) |
| 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 |
| 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 |
| 1239 |
|
|
| 1240 |
|
{\tt moveA:} |
| 1241 |
|
\begin{align*} |
| 1242 |
< |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 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) \\ |
| 1246 |
> |
\chi(t)\right), \\ |
| 1247 |
|
% |
| 1248 |
|
{\bf r}(t + h) &\leftarrow {\bf r}(t) |
| 1249 |
< |
+ h {\bf v}\left(t + h / 2 \right) \\ |
| 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) \\ |
| 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) \\ |
| 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) |
| 1261 |
> |
{T_{\mathrm{target}}} - 1 \right) . |
| 1262 |
|
\end{align*} |
| 1263 |
|
|
| 1264 |
|
Here $\mathrm{rotate}(h * {\bf j} |
| 1279 |
|
{\tt moveB:} |
| 1280 |
|
\begin{align*} |
| 1281 |
|
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1282 |
< |
\left\{{\bf j}(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) \\ |
| 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) \\ |
| 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) |
| 1296 |
> |
\chi(t + h) \right) . |
| 1297 |
|
\end{align*} |
| 1298 |
|
|
| 1299 |
|
Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to caclculate |
| 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 $h$ or $\tau_T$ can result in non-conservation |
| 1316 |
|
of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the |
| 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 |
|
|
| 1390 |
|
|
| 1391 |
|
{\tt moveA:} |
| 1392 |
|
\begin{align*} |
| 1393 |
< |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 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\} \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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\} \\ |
| 1420 |
> |
- {\bf R}_0 \right] \right\} ,\\ |
| 1421 |
|
% |
| 1422 |
|
\mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)} |
| 1423 |
< |
\mathsf{H}(t) |
| 1423 |
> |
\mathsf{H}(t). |
| 1424 |
|
\end{align*} |
| 1425 |
|
|
| 1426 |
|
Most of these equations are identical to their counterparts in the NVT |
| 1433 |
|
h / 2$. Reshaping the box uniformly also scales the volume of |
| 1434 |
|
the box by |
| 1435 |
|
\begin{equation} |
| 1436 |
< |
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)} |
| 1436 |
> |
\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}. |
| 1437 |
|
\mathcal{V}(t) |
| 1438 |
|
\end{equation} |
| 1439 |
|
|
| 1445 |
|
{\tt moveB:} |
| 1446 |
|
\begin{align*} |
| 1447 |
|
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1448 |
< |
\left\{{\bf j}(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\} \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) |
| 1469 |
> |
\chi(t + h) \right) . |
| 1470 |
|
\end{align*} |
| 1471 |
|
|
| 1472 |
|
Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required |
| 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 \cdot \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) \\ |
| 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} |
| 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 |
|
|
| 1537 |
|
|
| 1538 |
|
{\tt moveA:} |
| 1539 |
|
\begin{align*} |
| 1540 |
< |
T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\ |
| 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\} \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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\} \\ |
| 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)} |
| 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. |
| 1581 |
|
{\tt moveB:} |
| 1582 |
|
\begin{align*} |
| 1583 |
|
T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\}, |
| 1584 |
< |
\left\{{\bf j}(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\} \\ |
| 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) \\ |
| 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) \\ |
| 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) \\ |
| 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) |
| 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 |
| 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 |
| 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}} |
| 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}} |
| 1987 |
> |
\caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}} |
| 1988 |
|
\label{oopseTb:gofrs} |
| 1989 |
|
\begin{center} |
| 1990 |
|
\begin{tabular}{|l|c|c|} |
| 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 |
|
|
